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o 

We present the basic concepts and our recent developments in the density functional 
. approaches with the Skyrme functionals for describing nuclear dynamics at low energy. 

1 The time-dependent density-functional theory (TDDFT) is utilized for the exact lin- 

ear response with an external perturbation. For description of collective dynamics 
beyond the perturbative regime, we present a theory of a decoupled collective sub- 
manifold to describe for a slow motion based on the TDDFT. Selected applications are 
shown to demonstrate the quality of their performance and feasibility. Advantages and 
disadvantages in the numerical aspects are also discussed. 
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1. Introduction 

1.1. Nucleus as a quantum object 

The nuclei provide the matters with mass, the stars with fuel, and the universe with a 
variety of elements. It was discovered by Ernest Rutherford and coworkers about 100 years 
ago [1], that explains the large-angle scattering of alpha particles by a gold foil [2]. This 
discovery was also a beginning of the era of the quantum mechanics. Rutherford estimated 



the upper limit of the nuclear size which turned out to be much smaller than that of the 
atom. In the atomic scale (A), the nuclear size (fm) could be regarded as just a point! 
According to the classical mechanics, the atoms must collapse into nuclear size, because 
the attractive Coulomb potential eventually brings all the electrons into the nucleus. This 
mystery stimulated Niels Bohr to develop his idea on the quantum mechanics [3]. 

1.1.1. Atoms and molecules. With a knowledge of the quantum mechanics, it is easy to 
understand why the atoms do not collapse to the nuclear scale. If an electron was confined 
in the nuclear scale of femtometer, the uncertainty principle immediately tells us that its 
zero-point kinetic energy would become gigantic (order of GeV). In order to decrease this 
kinetic energy to a reasonable magnitude, the electron's wave function must have an atomic 
size of A. Thus, the atomic size is a consequence of the quantum effect. 
The molecules (and the solids) are a good contrast to atoms. The atom is bound by the 
Coulomb interaction, whose range is infinite (Vjv-e ~ ~Ze 2 /r), between a positively charged 
nucleus and electrons with a negative charge. Since the molecules consist of these charge- 
neutral atoms, the interaction between a pair of neutral atoms does not have the long-range 
tail of 1/r, but normally has a short-range repulsive part and an intermediate -rein £e cvttr ac- 
tive part (Fig. 1(a)). Because of this charge neutrality, the molecule is easy to disintegrate 
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Fig. 1 Schematic representation of atomic and nucleonic interactions V(r) as a function 
of the relative distance r. The energy of the bound (ground) state E is shown by a horizontal 
line, (a) Atom-atom interaction: Vo = 1 ~ 10 eV and c ~ as, where as is the Bohr radius, (b) 
Nucleon-nucleon interaction: V$ ~ 100 MeV and c = 0.5 ~ 1 fm. The essence of the present 
figure is taken from Fig. 2-36 in Ref. [4]. 

into smaller units (molecules and atoms). The atomic size is approximately constant and 
independent from the atomic number, while the molecular size varies depending on the num- 
ber of atoms and their kinds. Last but not the least, in the zero-th order approximation, 
the ground states of the molecules can be classically described as atoms located at fixed 
relative positions. The hindered quantum fluctuation in molecules is simply due to the fact 
that the atomic mass, which is approximately identical to the nuclear mass, is about 2,000 
times larger than the electronic mass. 

1.1.2. Nuclei. The nucleus has a number of properties analogous to the molecules, except 
for its strong quantum nature. It is a self-bound system composed of fermions of spin 1/2 and 
isospin 1/2 with approximately equal masses, called nucleons (protons and neutrons). The 
nuclear species are classified by the numbers of neutrons (N) and protons (Z). Rutherford 
discovered that the size of the nucleus is as tiny as the femtometer, but later it was found 
that the nuclear size varies, as its volume is roughly proportional to the mass number (A = 
N + Z). Each nucleon is a color-singlet (neutral) object. The interaction between a pair 
of nucleons (nuclear force) has a finite range of a scale of the pion Compton wave length 
(A,,-). Similar to the molecules, it has a short-range repulsive part and an intermediate-range 
attractive part (Fig. 1(b)). The nucleus can be disintegrated into small pieces with a small 
separation energy. In fact, heavy nuclei can have "negative" separation energies because of 
the repulsive Coulomb energy among protons. 

The quantum nature is an important difference between the nucleus and the molecule. In 
Fig. 1, we show schematic pictures of the atomic and nucleonic potentials. The ground state 
of the diatomic molecule (panel (a)) is formed at the bottom of the potential, E ~ — Vq. In 
the length scale of A, the atomic mass is heavy enough to localize the wave function at the 
location of the bottom of the potential, r « c. Thus, the relative distance between a pair of 
atoms is fixed at r = c. This property allows us to describe the atomic motion in the classical 
mechanics, such as in the molecular dynamics. In contrast, the nuclear interaction is not 
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strong enough to bind nucleons at the bottom of the potential. In other words, the nucleon's 
mass is too light to localize the wave function in the sub-femtometer scale. In deuteron, 
the zero-point kinetic energy cancels the negative potential energy ((T + V) « 0), leading 
to a bound state at approximately zero energy (E = —2.2 MeV 3> — Vo). The deuteron wave 
function spatially extends beyond the range of nuclear force (~ A,,-), which reduces the kinetic 
energy (T). This shows a striking contrast to the diatomic molecule, and is even analogous 
to the atoms, that the large size of the deuteron is a consequence of the quantum effect. This 
strong quantum nature also tells us that the infinite nuclear matter will not be crystallized 
even at zero temperature, but will stay as the liquid. The nuclear system provides us with 
unique opportunities to study femto-scale quantum liquids. 

1.2. Computing nucleus from scratch 

The strong quantum nature in finite nuclei leads to a rich variety of unique phenomena. 
Remarkable experimental progress in production and study of exotic nuclei requires us to 
construct theoretical and computational approaches with high accuracy and a reliable pre- 
dictive power. Extensive studies have been made for constructing theoretical models to 
elucidate basic nuclear dynamics behind a variety of nuclear phenomena. Simultaneously, 
significant efforts have been made in the microscopic foundation of those models. 
For light nuclei, the "first-principles" large-scale computation, starting from the bare 
nucleon-nucleon (two-body & three-body) forces, is becoming a current trend in theoret- 
ical nuclear physics. Among them, the Green's function Monte Carlo (GFMC) method is 
the most successful first-principles approach to nuclear structure calculation [5]. In this 
approach, using the Monte Carlo technique, the many-body wave function is sampled in the 
coordinate space with spin and isospin degrees of freedom. The success of the GFMC clearly 
demonstrates that we are able to construct a light nucleus from the scratch on the computer. 
The GFMC method has been applied to nuclei up to the mass number A « 10. Another first- 
principles approach is to project the nuclear Hamiltonian in a truncated Hilbert space, then 
diagonalize it. This is called no-core shell-model (NCSM) method [6]. The NCSM also shows 
successful applications up to the p-shell nuclei. The GFMC and NCSM both indicate the 
exponential increase in computational tasks with respect to the increasing nucleon number. 
The third approach, the coupled-cluster method (CCM), has an advantage that the required 
computation increases only in power with respect to the nucleon number. The CCM, which 
was originally invented in nuclear physics [7] and later became extremely successful in quan- 
tum chemistry [8], has been revisited as an ab-initio computational approach to nuclear 
structure [9]. Especially, it is powerful to study closed-shell nuclei. 

Although these first-principles approaches have recently shown a significant progress, they 
are still limited to nuclei with the small mass number. This may sound mysterious to physi- 
cists in other fields, because we know that the similar kinds of approaches are able to 
treat systems of much larger particle numbers. For instance, using the CCM, nowadays, 
the chemists can easily calculate a molecular structure with 100 electrons. Why is the first- 
principles calculation of nuclear structure so difficult? The answer is perhaps trivial for 
nuclear physicists, but may not be so for others. Let us pick up several important aspects 
leading to this answer. (1) The strong quantum nature. As we have discussed previously, 
the full quantum mechanical treatment is necessary for nuclear structure calculation. (2) 
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Strong coupling nature. The nucleon-nucleon scattering length is approximately a « — 18 
fm in 1 So channel. This is much larger than the mean distance between nucleons inside the 
nucleus (kp\a\ S> 1). (3) Singular property of the nucleonic interaction. To reproduce the 
phase shift for the nucleon-nucleon scattering, the interaction should have a strong repulsive 
core at short distance [4]. (4) Complexity of the nucleonic interaction. The interaction has 
a strong state dependence which may be represented by the spin- and isospin-dependence 
[4]. It also contains strong non-central and non-local terms, such as tensor and spin-orbit 
interactions. Furthermore, for a quantitative description of nuclei, it is indispensable to intro- 
duce three-body interactions in addition to the two-body force. (5) Coexistence of different 
interactions. In addition to the strong interaction, we need to treat the electromagnetic 
interactions among protons. (6) Coexistence of different energy scales. The nuclear binding 
energy amounts to order of GeV for heavy nuclei. Because of the strong quantum nature 
in nuclei, this is a consequence of the cancellation between positive kinetic energy and neg- 
ative potential energy. Thus, we need to compute these enormous positive and negative 
components in high accuracy to obtain a correct binding energy. (7) Fermionic nature of 
nucleons. Needless to say, since the nucleons are fermions, the total wave function must be 
anti-symmetrized. (8) Finite systems without an external potential. Electronic many-body 
problems in molecules and solids are solved with external fields produced by nuclei with 
positive charges. In contrast, the nucleus is a self-bound finite system. Thus, we are usually 
required to obtain an intrinsic wave function without the center-of-mass degrees of freedom. 
This requirement often restricts our choice of the basis functions. 

Despite of these difficulties, significant advances in the computer power may lead to a realistic 
"first-principles" construction of the sti-shell nuclei in near future. There have been extensive 
efforts toward this direction, especially developing an algorithm suitable for parallel use of 
the vast number of processors [10]. 

1.3. Density functional theory (DFT) 

In contrast to the first-principles calculations, which are limited to nuclei with small number 
of nucleons, the density functional theory (DFT) is currently a leading theory for describing 
nuclear properties of heavy nuclei [15, 16]. It is capable of describing almost all nuclei, includ- 
ing nuclear matter, with a single universal energy density functional (EDF). An argument 
based on the quantum many-body theory leading to nuclear EDF was also developed in 
70's— 80's [11, 12]. In addition, its strict theoretical foundation is given by the basic theorem 
of the DFT [13, 14]. Since the nucleus is a self-bound system without an external potential, 
we should slightly modify the DFT theorem. This will be discussed in Sec. 2. 
The nucleus by itself produces a potential confining nucleons which is analogous to the Kohn- 
Sham potential in DFT [14]. Nuclear physicists often call this potential "mean field", though 
it is different from a naive mean-field potential directly constructed from the nucleonic 
interaction. There are many evidences for the fact that the mean-free path of nucleons inside 
the nucleus is significantly larger than the nuclear radius [4], in spite of the strong (and even 
singular) two-body interaction. This is partially due to the Pauli exclusion principle. Since 
the nucleon mean-free path is roughly proportional to {ei?/(e — £f)} 2 [17], it is significantly 
enhanced for nucleons whose energies are close to the Fermi energy. Another even simpler 
argument was given by Bohr and Mottelson [4], that the nuclear normal density is much 
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lower than the one giving the close packing limit (crystalline limit). In this argument, the 
quantum effect plays a primary role. 

The DFT theorem guarantees the existence of generalized density functionals for every 
physical observable (see Sec. 2). However, to construct the exact functional, we need exper- 
imental data and other theoretical inputs. Currently, there are a variety of EDFs which 
predict somewhat different properties for nuclei very far from the stability line. We are in 
want of finding an ultimate universal energy density functional, which is capable of exact 
description of every nucleus in the nuclear chart. In addition to the recent progress in the 
first-principles calculations for light nuclei, the radioactive beam facilities in the world will 
give us ideal opportunities to determine the parameter set for a better functional. New data 
on neutron halos and skins in medium heavy nuclei may provide important information 
on its dependence on density and density gradient. Observation of new isotopes in a long 
isotopic (isotonic) chain may lead to useful constraints on the isovector parts of the energy 
functional. 

An extension of the DFT to the time-dependent DFT (TDDFT) provides a feasible descrip- 
tion of many-body dynamics, which contains information on excited states in addition to the 
ground state. The TDDFT is justified by the one-to-one correspondence between the time- 
dependent density and time-dependent external potential [18], which will be presented in 
Sec. 2.4. The TDDFT has vast applications to quantum phenomena in many-body systems. 
Among them, the perturbative regime has been mostly studied so far. Different approaches 
to the linear response calculations will be presented in Sec. 3. It is of significant interest but 
challenging to go beyond the perturbative regime. Nuclei show numerous phenomena related 
to the large amplitude collective motion, such as fission, shape transition, shape coexistence, 
anharmonic vibrations, and so on. We present, in Sec. 4, a theory to identify an optimal 
collective submanifold in the classical phase space of large dimensions. 

2. Basic formalism for particles in self-bound systems 

The density functional theory (DFT) has been extremely successful for calculations of 
ground-state properties of atoms, molecules, and solids. It describes a many-particle system 
exactly in terms of its one-body density alone. The DFT is based on the original theorem of 
Hohenberg and Kohn (HK) [13] which was proved for the ground-state of the many-particle 
system. Every observable can be written, in principle, as a functional of density. 
In nuclear physics, however, we need to treat an isolated system without an external poten- 
tial. The present nuclear EDF produces a localized density profile without an external 
potential. Thus, the ground state spontaneously violates the translational symmetry and 
seems to contain spurious excitation related to the center-of-mass motion. Furthermore, it 
also violates the rotational symmetry when the ground state is spontaneously deformed. It 
is interesting to see whether the nuclear EDF can be theoretically justified in a strict sense. 
We present a possible justification in Sec. 2.1, according to the recent progress [20-22]. 

2.1. DFT theorem for a wave-packet state 

The HK theorem [13] guarantees one-to-one mapping between a one-body density p(r) and 
an external potential vo(r). Then, since the ground-state wave function is a functional of 
density, in principle, all the observables should be functionals of the density as well. However, 
to describe an isolated self-bound finite system in a box of volume V, it is somewhat useless to 
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use the ground-state density in the laboratory frame, because it must be constant, p(r) = p = 
N/V — > (V — > oo), where N is the particle number. Instead, we want to use a functional of 
the intrinsic density p(f — R), where R is the center-of-mass coordinate of the total system. 
In this case, the original HK theorem cannot be justified, because it adopts a one-body 
external field vo(f) coupled to the density p{r) in the laboratory frame. 
Validity of the DFT for the intrinsic states has been discussed by several authors [20-22]. In 
this section, we present a method to define the DFT for an intrinsic state, more precisely for 
a "wave-packet" state. The argument here essentially follows the idea by Giraud et al [22]. 
In fact, all the nuclear EDFs, currently available, produce a wave-packet state. The mini- 
mization of an EDF E[p] without the external potential (vo(f) — > 0) leads to the nucleon 
density distribution p(r) with a finite radius. This violates the translational symmetry of 
the original Hamiltonian. Therefore, it is of our interest here to justify the energy functional 
of the wave-packet density in the laboratory frame in a strict sense. 

First, we assume that a wave-packet state in the laboratory frame can be expressed by a 
product wave function of intrinsic and spurious degrees of freedom, |$) = \(f>) (g) |%). Here, 
\(p) indicates the intrinsic state and |x) defines the spurious motion. This decomposition can 
be exactly done for the translational motion. 

n = H^,«) + J^, (i) 
Hn,--- ,f N ) = ,6v_i)®xCR), (2) 

where R (P) denote the center-of-mass coordinates (total momenta). & and 7Tj are the 
relative Jacobi coordinates and their conjugate momenta, respectively. Since the intrinsic 
ground state, which is supposed to be unique, is completely independent from the spurious 
motion, we can adopt an arbitrary form of x(-R); e.g., x(R) exp(— R 2 /2b 2 ). Then, the 
ground wave-packet state can be obtained by the variation after the projection: 
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where the projection operator P does not change the intrinsic state \(f>) but makes x(R) an 
eigenstate of the total linear momentum P = 0. The variation with respect to the full space 
(|$) = \4>) <8> |x)) contradicts the uniqueness of the ground state, because states with different 
x(R) give the same P\&). Therefore, the variation here should be performed with respect 
only to the intrinsic state \<fi). This is indicated by the subscript "x:fixed" in Eq. (3). The 
wave-packet density profile is simply given by p(r) = (^^(r)^) ; that depends on the 
form of x(R)- I n the followings, we always assume a fixed form of x(R) f° r t ne wave-packet 
state |$). 

Next, we introduce an external potential Vo = X^ofji)- The following minimization with 
respect to the intrinsic state \4>) leads to the "minimum" energy Ebfao] an d defines the wave 
packet |$). 
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Noted that Vq operates on a state |3>), not on a projected state P\&). Eq[vq] does not 
correspond to the ground-state energy of a system with the Hamiltonian H + Vq, however, 
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it reduces to Eq. (3) for Vo = 0. Now, let us show the one-to-one correspondence between 
the external potential Vq and the wave-packet density p(r). The proof proceeds by reductio 
ad absurdum in the same manner as the original proof of HK [13]. Assume that another 
potential v' (f), which defines the wave packet |3>'), produces the same density p(f). Then, 
the energy for |<J>') is given by 

(ti>'\HP\&') 

^o] = { { ^ ) ) + (^Km- (5) 

If we replace the state |$') by |<E>), the energy must increase. 
Eo[v' ] < 

($\HP\$) 

= -jk\PW + {m ^ ) + mv "~ v ° m 

= E [vo] + J p{r){v' Q {r)-v Q {r)}d z r. (6) 

Interchanging Vq and Vq, we also find 

Eq[v ] < Eq[v' ] + J p(r) {v (r) - v' (r)} d 3 r. (7) 

Addition of Eqs. (6) and (7) leads to the inconsistency, i^o^o] + -^oK)] < Eq[vq] + Eq[v'q}. 
This proves the one-to-one correspondence between the external field Vq and the wave-packet 
density p{r). Thus, both vo(r) and the wave packet |$) are functionals of the density p(r). 
In order to lift restriction to ^-representative densities, we use the constrained search [23] in 
which one considers only states that produce a given density p(r), and define the universal 
functional 

'{&\HP\&) 



F x\p\ = ™ in 

<3>— >p 



(8) 



<$|P|$) . ' 

Here, the subscript "<£ — > p" indicates the minimization with a constraint on p(r). The 
density functional F x [p] contains the energy of the fixed center-of-mass spurious motion, 
E x , that is trivially given by E X (N) = {x\P 2 \x)/(?Nm). Therefore, the energy of the ground 
state with P = (intrinsic energy) may be obtained by the minimization with a constraint 
on the total particle number, as 



E- mtr (N) = E p= q{N) = min 

p 



F x [p]-\^ f p(r)d 3 r-N^ 



E X {N). (9) 



Note that, although the wave function |%) is fixed, E x may depend on the total mass of the 
system (particle number N). In principle, the expectation value of any observable O, which 
only depends on the intrinsic degrees of freedom, is a functional of p, 

o[ P ] = (01010) x ( x \x) = mom, (io) 

because the wave-packet state |<&) is given as a functional of p. Since the state |%) is fixed, 
there is a trivial correspondence between the wave-packet density p(r) and the intrinsic 
density p(r — R). This completes the basic theorem of the DFT for the wave packet. 
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2.2. The Kohn-Sham (KS) scheme 

For a many-body system of fermions, the shell effects play a major role to determine the 
ground state. In other words, we need a density functional which takes account of the kinetic 
energy properly. This is known to be difficult in the local density approximation [19]. At 
present, a scheme given by Kohn and Sham [14] only provides a practical solution for this 
problem. Here, we follow the same argument. 

We introduce a reference system which is a "virtual" non-interacting system with an external 
potential v s {r). This reference system is supposed to reproduce the same density p(r) of the 
"physical" interacting wave packet, but does not have to reproduce the center-of-mass wave 
function |%). The ground state of the reference system is trivially obtained as a Slater 
determinant constructed by the solution of 1 

-^V 2 + v s (r^j fa{r) = ei^i(r), (H) 

adopting the unit h~ = 1, and the density is given by p{r) = Zti \Mr)\ 2 - The kinetic energy 
is given by 

N ( 1 \ 

r.M = XW-^v 2 J|*>. (12) 

i=i ^ ' 

The variation of the total energy of the reference system 

E s [p] = T s [p] + f v s (f)p(r)d 3 r, (13) 

with a constraint on the particle number, S(E s [p] — p f p(r)d 3 r) = 0, leads to the following 
equation: 

Although we did not explicitly construct T s [p] as a functional of p(r), the solution of Eq. (14) 
must be identical to the solution of Eqs. (11) and (12). 

The success of the Kohn-Sham (KS) scheme comes from a simple idea to decompose the 
kinetic energy in the physical interacting system into two parts; T s [p], which is a major 
origin of the shell effects, and the rest, which is treated as a part of "correlation energy" 
described by a simple functional of density, 

F x [p]=T s [p} + E c [p], (15) 

where E c [p] = F x [p] — T s [p\. Then, the variation of F x [p] leads to Eq. (14) but the potential 
is now a functional of density, defined by v s (f) = 5E c [p]/5p(r). The only practical difference 
between the reference system and the interacting system is that, since v s (r) is a functional 
of density in the latter, Eq. (11) must be self-consistently solved. These equations are called 
Kohn-Sham (KS) equations. The self-consistent solution of the KS equations provides the 
density p(f) of a wave-packet state with a fixed \\) corresponding to a (local) minimum of 
the EDF, F x [p]. The success of the KS scheme is attributed to the goodness of the local 
density approximation for E c [p\. 



1 Precisely speaking, the orbitals fa are not necessarily the eigensolutions of Eq. (11), but arbitrary 
as far as they give the same Slater determinant. We come back to this gauge freedom in Sec. 2.5. 

2 The HK theorem guarantees that T s [p] of the reference system is a functional of the density. 
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2.3. Open issues 

2.3.1. Subtraction of the center- of-mass energy. In the proof of the basic theorem for the 
wave packet, given in Sec. 2.1, we need to fix a wave function of the center-of-mass motion 
|x). The choice of this spurious wave function is arbitrary, but the energy E x depends on 
this choice. In practice, the subtraction of E x is normally performed by constructing the 
state |x) from the obtained result. This is somewhat inconsistent with the assumption of the 
fixed center-of-mass state |x). This could be easily corrected by taking E x of a given 
However, as far as we know, this has not been examined yet. 

2.3.2. Validity of the Kohn-Sham scheme. The KS scheme is to take into account a major 
part of the kinetic energy as T s [p], and the rest as a correction. In other words, the KS scheme 
implicitly assumes that the energy E c [p] is able to be well approximated by a simple func- 
tional of p. In fact, this question is still an open issue, not only in the nuclear physics but 
also in other quantum many-body systems. In the present wave-packet theory, the kinetic 
energy of the wave packet, T x [p] = (3>|T|<J>), depends on the center-of-mass state x(R)- There- 
fore, there may be an optimal choice for x(i?) to minimize the difference T x [p] — T s [p\. The 
question about the validity of the Kohn-Sham scheme remains to be answered. 

2.3.3. Non-spherical wave packet. The nuclear EDFs are known to produce a spontaneous 
symmetry breaking about the rotational symmetry. Namely, we often encounter a deformed 
wave-packet density, which accounts for appearance of the rotational spectra in nuclei. For 
instance, many experimental evidences indicate that nuclei in the rare-earth region and in 
the actinide region are deformed [24]. According to the argument in Sec. 2.1, we may separate 
the rotational motion from the intrinsic degrees of freedom, then, we have 

$(ri,--- ,nvW(£)®x(M), (16) 

where indicates angle variables. Then, replacing the operator P by that of angular momen- 
tum projection, the DFT for deformed wave packet can be shown in the same manner. 
However, in contrast to the translational motion, the separation of the rotation degrees 
of freedom is not exact, but only approximate. Thus, there remains an ambiguity for the 
definition of the functional (8): Namely, the minimization must be performed in the entire 
space except for the subspace that accounts for translational and rotational correlations. 
In this sense, the use of the EDF, which produces a deformed state, can be justified only 
approximately. 

2.4- Time- dependent density functional theory 

The DFT is designed for calculating the ground-state properties. For excited-state properties 
and reactions, the time-dependent density functional theory (TDDFT) is a powerful and 
useful tool. In this section, we recapitulate the basic theorem for the time-dependent density 
functional theory (TDDFT). 

Since the proof of the HK theorem is based on the Rayleigh-Ritz variational principle, its 
extension to the time-dependent density is not straightforward. This was done by Runge 
and Gross [18], showing that there is one-to-one correspondence between a time-dependent 
density p(r,t) and a time-dependent external potential v(f,t). The external potential is 
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required to be expandable in a Taylor series about the initial time to, 



k=0 



(17) 



The external potentials, v(f,t) and v'(r,t), are defined to be different if there exist some 
minimal nonnegative integer k such that Vu>fc(r) / where u>fc(r) = Vk{r) — v' k {r). In other 
words, v (r, t) and v'(f, t) differ more than a time-dependent function, v(r, i) — v'(f, t) / c(i). 
Note that the potentials differing by the time-dependent constant c(t) produce the same 
density p(f, t) because the corresponding wave functions differ by a merely time-dependent 
phase, as in Eq. (27) with a(t) = /* c(s)ds. 

Now, let us assume that starting a common initial state |^o) = |^(*o))) two different external 
potentials, v(f,t) and v'(f,t), produce densities p(f,t) and p'(f,t), respectively. From this, 
we first prove that the current densities, j(f,t) and j'(f,t), are different. Using the current 
density operator j(r), the equation of motion is written as 



ip(r,t) = ^(t)\[j(r),H(t)]\^(t)}, 



(18) 



where H{t) = Ho + X)£Li v ifh We have the same equation for j'(f,t), with |\P(t)) and 
v(f,t) replaced by | and v'(f,t), respectively. Then, we have 



d_ 

di 



{l(r,t)-f(f,t)} 



t=to 



- -i{V \\)(rl,H(to) - H'(to)}\Vo) 



= p(r,to)Vw (r). 

m 



(19) 



If Vwo{r) 7^ 0, it is easy to see that j(r,t) and j'(f,t) are different at t > to- In case that 
Vwo(f) = and Vwi(f) / 0, we need to further calculate derivative of Eq. (18) with respect 
to t. 

.d 2 



t=t 



(*o|b>1,^^]|*o) + (*o|[[J(r1,^(to)],^(to)]|*o). (20) 



The second term of Eq. (20) vanishes for d 2 / dt 2 {j (r, t) — j'(r,t)}\t=t , because H'(to) 
H(t ) + const. Thus, 



dt 2 



{l(r,t)-f(r,t)} 



t=t 



—p(f,to)V Wl (r)^0. 
m 



(21) 



Again, we can conclude that j(f,t) / j'(f, t) at t > to- In general, if Vwk(r) = for k < n 
and Vw n {f) / 0, we repeat the same argument to reach 



n+l 



}(f,t)-j'{r,t)} 



t=t 



= --p(f,t )Vw n {r)^0. 
m 



(22) 



Therefore, there exists a mapping from the expandable potential v(f, t) to the current density 
j(r,t). 
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Next, we use the continuity equation 

| {p(f, t) - p'(f, t)} = - V • {j(r, t) - j'(f, t)} , (23) 
and calculate the (n+l)-th derivative of Eq. (23) at t = to- From Eq. (22), 

n+2 



-Q t ) {p(r,t)-p'(f,t)} 



= -V-{p(f 1 t )Vw n (r)}. (24) 
m 

t=to 



Provided Vw n (r) / 0, we can prove that the right-hand side of Eq. (24) does not vanish 
identically. This is done by using the following identity: 

V • {p{f, t )w n (r)Vw n (r)} - w n {r)V • {p{f, t )Vw n (r)} = p{r, t ) {Vw n {r)} 2 . (25) 

The integral of both sides of Eq. (25) over the entire space leads to 

- J dfw n (r)V ■ {p(r, to)Vw n (f)} = J drp(f, t ) {Vw n (r)} 2 > 0, (26) 

where we assume p(f, to) is localized in space so that the surface integral vanishes. Therefore, 
from Eq. (24), we can conclude that the densities p(r,t) and p'(f,t) are different at t > to- 
This completes the proof. 

The time-dependent density determines the time-dependent external potential except for 
the time-dependent constant. Therefore, the many-body time-dependent state should be a 
functional of density except for a time-dependent phase. 

|*(t))=exp(-ia(t))|*[p](t)). (27) 

Any observable quantity must be independent from the global phase, a(i), thus, a unique 
functional of density, yl[p(t)]. Note that these time-dependent density functionals depend on 
the initial many-body state |^o)- 

2.5. Time- dependent Kohn- Sham (TDKS) equations 

In practice, we use the Kohn-Sham scheme [14] for numerical calculations. Assuming the 
u-represent ability, the time-dependent Kohn-Sham (TDKS) equation is given by 

i^i{r,t) = \-^-V 2 + v s [ P }(f,t))^(r,t), i = !,-••, N. (28) 



The density of a system is expressed by p(f, t) = YliLi iV'j^ *)| 2 - I n practice, we usu- 
ally adopt the potential same as the one for calculation of the ground state ("adiabatic 
approximation"), except for the external potential v(f, t). 



v s [p](r,t) = v{r,t)+ 6Lj>] 



5p(f) 



(29) 



The density is invariant with respect to the unitary transformation U (N) among occupied 
KS orbitals. Therefore, there are gauge degrees of freedom to choose this transformation 
at any instant of time. For explicit notification of the gauge freedom, it is convenient to 
introduce the matrix notation as follows. Let {\a)} be an arbitrary single-particle basis set, 
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and we define a matrix ip of size of oo x N, as tp a i(t) = (<x\tpi{t)}. Then, the density matrix 
p(t) can be written as 

PaP {t) = j2(<*Mt))mw) = (mftw) R ■ (30) 

i 

The orthonormal property of the KS orbitals is expressed as ^ (t)^{t) = 1. Denoting the 
TDKS Hamiltonian in Eq. (28) as h s (t) = h s [p{t)\ = -V 2 /(2m) + v s [p](f, t), the TDKS 
equations (28) can be generalized into the following form. 

i^m = hsW(t)-i>(t)Z(t), (31) 

where £(t) is an arbitrary N x N time-dependent Hermitian matrix which represents a gen- 
erator of the U (N) transformation. Equation (31) is equivalent to the well-known equation 
for the density matrix [25]. 

d 

i-Q t P(t) = [h s (t),p(t)}. (32) 
The stationary state corresponds to the time-indenpendent density, dp/dt = 0. 

2.6. Pairing correlations: Kohn-Sham-Bogoliubov (KSB) equations 

The HK theorem (or its wave-packet version), in principle, guarantees that the energy of 
the system can be exactly written as a functional of density, E[p\. However, in practice, it 
is often difficult to take into account all the correlation energy in E[p], solely depending 
on p(r). The kinetic energy is such an example, which is resolved by the genius idea by 
Kohn and Sham. The pairing correlation energy E pa [ T , which is important for heavy nuclei 
in open-shell configurations, is another example difficult to be expressed by p(f) only. 
To overcome this difficulty, a common strategy is to extend the KS equations, according to 

the Bogoliubov's quasiparticles [26]. Each orbital now has two components, Q v = i^y^j ' an< ^ 

its number is basically infinite (y = 1, ■ ■ ■ , oo). These are called quasiparticle (qp) orbitals. 
The previous KS equations are extended to the following equations, which we call Kohn- 
Sham-Bogoliubov (KSB) equations 3 hereafter: 

{Us - MO *u = E u $„, (33) 

where 

«>-=(xu-%^ "-(;-)■ 

The Hamiltonian h is in the same form as that in Eq. (11), h = — V 2 /(2m) + v s [p, n], while 
the pair potential A[p, n] is introduced to describe the pairing correlations. The same form 
of equation is known as the Hartree-Fock-Bogoliubov equation in nuclear physics [19, 27]. 
Now, the KSB Hamiltonian in Eq. (34) not only depends on the density p, but also on the 
pair density k. In Eq. (33), the matrix convention is assumed. Namely, when we adopt a 
single-particle basis of {|a)}, we have U au = {a\U u ) and V av = (a\V u ), and h a p and 
correspond to Hermitian and anti-symmetric matrices, respectively. 



; Again, this corresponds to a special choice in the gauge degrees of freedom. See Sec. 2.7. 
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The chemical potential u is determined so as to satisfy the following number condition: 
tip = N. With this number constraint, trp must have a finite value. On the other hand, 
the values of the pair density k are determined solely by the variation of the total energy. 
Therefore, the self-consistent solution of the KSB equations (33) spontaneously produces the 
finite values of k and A. 

The solutions of the KSB equation have a "paired" property: If the qp state <& v = ( is 

( y*\ 

a solution of Eq. (33) with eigenvalue E u , the qp state <& v = I ^ I is also a solution with 

eigenvalue —E v . We call $^ "unoccupied" qp orbitals, and Q v "occupied" qp orbitals [28]. 
This naming is based on the generalized density matrix, 



*= ' , -. <»> 




which is Hermitian and idempotent: R 2 = R. The "unoccupied" ("occupied") qp orbitals 
correspond to the eigenvectors of R with eigenvalue (1); R<& v = and R<h v = <& v . 
Denoting the dimension of the single-particle Hilbert space as M, we may define the 2M x M 
matrix <3? as follows. 

= ( a ^-- U (30) 

{{a - M\V V ) a = M + !,■■■, 2M 

which represents "unoccupied" qp orbitals (y = 1, ■ ■ ■ , M). The "occupied" orbitals <l with 
size of 2M x M are defined in the same manner, with U u (V v ) replaced by V* (U*). The 
generalized densities are expressed in terms of these matrices as R = = 1 — . The 
orthonormal property of the qp orbitals is given by = = 1. Combining the 'unoc- 
cupied" and "occupied" orbitals to construct the 2M x 2M matrix W = (<£,$), the matrix 
W becomes a unitary matrix [19]. 

In the energy functional of Skyrme type [29], the pairing correlation energy is simply added 
to the original energy functional. 

E[ P ,K} = E[p] + E pail [p,Ki (37) 

that depends only on the local densities. Therefore, Eq. (33) becomes local in coordinates. 
However, in general, the KSB Hamiltonian, h and A, are not necessarily local. For instance, 
the Gogny functional [30] gives non-local KSB equations. 

2.7. Time- dependent Kohn-Sham-Bogoliubov (TDKSB) equations 
For a time-dependent description, the inclusion of the pair density leads to the time- 
dependent Kohn-Sham-Bogoliubov (TDKSB) equations. They can be formulated in a matrix 
form analogous to Eq. (31). Using an arbitrary 2M x M Hermitian matrix E(i), the TDKSB 
equations may be written as 

i^(i) = W,(t)^)-«(()S(i), (38) 

where the TDKSB Hamiltonian is given by Eq. (34), which depends on time through the 
densities p(t) and Here, ^(t) represent time-dependent "unoccupied" qp orbitals [y = 
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1, ■ ■ ■ , M). The "occupied" orbitals ^(t) are defined in the same manner, with U u (V u ) 
replaced by V* (17*). The TDKSB equation (38) holds for *(t), as well. 
Analogous to the stationary case, the generalized density R(t) can be written as R(t) = 
^S>(t)^>^(t) = 1 — ty(t)^(t), and the orthonormal property of the qp orbitals is given by 
the unitarity of the 2M x 2M matrix W(i). The "unoccupied" ("occupied" *(t)) 

correspond to the subspace with eigenvalue (1); R(t)^f(t) = and R(t)^(t) = ^(t). In 
the generalized density matrix formalism, the TDKSB equation is written in an analogous 
manner to Eq. (32): 

i^R{t) = [H,(t),R{t)]. (39) 

So far, we have shown similarities between Eqs. (31) — (32) and (38) — (39). However, there 
is an important difference between Eq. (32) and Eq. (39). The stationary solution in Eq. 
(32) corresponds to dp/dt = 0. In contrast, it is not the case in Eq. (39), dR/dt ^ 0. Let us 
examine this difference in details. The TDKSB equation (38) can be recast into another form, 
convenient for taking its stationary limit. First, let us factor out the time-dependent phases 
as follows: *(£) = exp(— ipJ\ft)W(t) and *(£) = ex.p(—ipj\ft) *'(£). Here and hereafter, we 
denote the remaining parts of the quantities as the "primed" ones. The generalized density 
becomes 

R(t) = *(t)*t( t ) = exp(-ipNt)R'(t) exp(+ipj\ft), (40) 

where 

*w = Wn = ($l !_%)' (41) 

Namely, the transformation *(i) — > *'(£) does not change the density p, but modifies the pair 
density as K,(t) = e~ 2tflt k' (t) . Since the pair potential A(i) is usually a linear functional of 
K(t), the same time-dependent phase should be assumed for A(i) as well: A(t) = A'(i)e _2v *. 
The Hamiltonian is transformed in the same manner: 



H a (t) = exp(-ipNt)n' s (t) exp(+ipMt). (42) 
With these primed quantities, the TDKSB equation (38) can be rewritten as 

i^'(t) = [H'M - pM) V(i) - V{t)E{t), (43) 

or equivalently, in the generalized density matrix, 

i^ t R'{t)=[U' s {t)-pN,R\t)]. (44) 

It is now clear that the stationary solution corresponds to dR'/dt = 0, not to dR/dt = 0, 
with a proper choice for the parameter p identical to the chemical potential. In Eq. (43), 
it corresponds to d^f'/dt = with a choice of the M x M gauge matrix E UI/ i = {^'^T-L'g — 
pN^'y,). It should be noted that the generalized density R{t) is invariant with respect to 
the choice of the gauge matrix In contrast, the time-dependent phase factor in re(i) 

and A(t) have a physical origin and cannot be removed by the gauge choice. In fact, it 
is a boost transformation, e - *^' , from the laboratory frame of reference to the body-fixed 
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frame. The stationary solution with k 7^ (A 7^ 0) corresponds to a time-dependent solution 
in the laboratory frame: 

( p -int n \ 
^ ,-•<••) (45) 

This is a collective motion associated with the spontaneous generation of the pair density, 
called pair rotation. Therefore, in terms of the TDKSB formalism, the appearance of the 
chemical potential in the stationary KSB equation (33) comes from the boost transformation 
to the body-fixed frame rotating in the gauge space. This is analogous to the appearance of 
the cranking term —ujJ x in the spatially rotating frame of reference [19]. In the case of pair 
rotation, since the particle number is finite N > 0, the system is rotating in the gauge space, 
even at the ground state. This rotation affects the intrinsic modes of excitation, thus, the 
Hamiltonian in the rotating frame, % s — fiAf, should be utilized to calculate the intrinsic 
excitation spectra. This point will be discussed again in Sec. 4.5. 

3. Perturbative regime: Linear response 

The theorem of the TDDFT tells us that the functional may depend on the initial state, 
I* (to)). This additional ambiguity can be removed by assuming that the initial state is 
identical to the ground state. With this assumption, the linear response theory with a weak 
time-dependent perturbation is formulated in this section. The formulation is basically iden- 
tical to the one known as the random-phase approximation in nuclear physics [19, 27]. 
However, according to the concept of the TDDFT, the theory gives the exact linear density 
response, with no approximation involved, in principle. 4 

3.1. Time- dependent linear density response 

We consider a system subject to a time-dependent external potential 

v(f,t) = { (46) 
[vi(r,t) t>0 

in addition to the static potential vo(f) of the unperturbed system. 5 In this section, we 
use the notation of the four vector x = (r,t). We assume that the system is at the ground 
state at times t < 0. Thus, the initial density po(r) at t < can be obtained from the self- 
consistent solution of the Kohn-Sham equations (11). The first-order density response, p{x) ~ 
Po(r) + Pi(x), is given by 

«(*) = /^<n(^>M (47) 



with the density-density response function 

8p(x) 



U(x,x') 



(48) 

v=0 



5v(x' 

The right-hand side of Eq. (48) is a well-defined quantity, since the basic theorem of TDDFT 
in Sec. 2.4 guarantees that the time-dependent density is a functional of the time-dependent 
external potential; p[v](a;). 



4 In practice, some approximations are involved, such as the adiabatic approximation of Eq. (29). 

5 For isolated nuclear systems, we have v — 0. 
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For non-interacting particles moving in an external potential of v s (x), there is a one-to-one 
correspondence between the time-dependent density and the potential. Therefore, we have 

p(x) = p[v s ](x), v s (x) = v s [p](x). (49) 

The density-density response function for the non-interacting system is given by 

5p(x) 



U s (x,x') 



8v s (x') 



Vs [po] 



(50) 



The potential v s (x) is written as a sum of a given external potential and the rest of the Kohn- 
Sham potential, v s {x) = v(x) + Uk s [/>](x). For instance, in the adiabatic approximation of Eq. 
(29), vy iS (x) = SE[p]/5p(x). Therefore, using the chain rules, Eq. (48) can be connected to 
its non-interacting U s (x,x'): 



U(x,x') 



d A y 



5p(x) 



o~v s (y) 



= J d 4 yU s (x,y) L(y - x') + J d A y § ^ 
= U s (x,x') + J d A y J d A y'U s (x,y)w(y,y')U(y',x'), 



W) 

po 5V ( X ') 



v=0 



where the residual kernel is given by 



w (x, x ) = 



Sv ks (x) 



(51) 



(52) 



Po 



5p(x') 

In the adiabatic approximation, this is equal to the second derivative of the energy functional 



w(x, x') = 



(53) 



5p{x)5p{x') 

Most of the energy functionals currently available are local in time, which leads to w(x, x') cx 
5{t-t'). 

Multiplying the Dyson-type equation (51) by the perturbing external potential v\(x) leads 
to the linear density response of Eq. (47). 

pi(x) = J d 4 x'U s (x,x')v sci (x'), (54) 

where the self-consistent effective field, given by 

consists of the external perturbation v\ and the induced residual field v^\{x) = 
J d 4 yw(x,y)pi(y). Thus, the exact representation of the linear density response pi(x) of 
a real interacting system can be written as the linear density response of a non-interacting 
system to the self-consistent effective perturbation v sc f(x). 

The formal solution for the density response II is given by solving the Dyson-type equation 
(51), II = (1 — II S ■ u>) _1 • II S . The non-interacting response II S is explicitly given in the 
followings, and the residual kernel w(x, x') is usually calculated using the adiabatic approx- 
imation of Eq. (53). In this response function formalism, the II (a; ) is usually solved in the 
frequency domain, to calculate the strength function, transition density, etc. 
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3.2. Linear density response with the Green's function method 
The Fourier transform brings Eq. (54) into 

Pi(r» = j d 3 r'U s (f,r';uj)v sci (r',uj), (56) 

where the frequency-dependent effective field is given by 

v sc{ (f,uj) = vi(f,u) + J d 3 r' w(f,r')pi(r',uj). (57) 

The non-interacting response function II S is expressed in terms of the static Kohn-Sham 
orbitals 4>k and their eigenenergies e&: 

tj ,-. -) \ / -a <Pm{r)4>* m (r') 

n s (r,r» = ^^{r)Ur') ^ e . +UJ _ e ■ ~ 

i<iv m>N 

(58) 

The restriction for the summation with respect to the index m can be lifted, because the 
first and second terms in Eq. (58) give the same magnitude but with an opposite sign for 
m < N. Using the spectral representation of the one-particle retarded Green's function for 
non-interacting particles, 

_£M2fi£>, (59) 
^ u - e k + in 

one may replace summed orbitals with respect to m in Eq. (58) by the Green's function. 

n 8 (f,r';w) = Y J {<t>l(^<t>i{r')G {+ X?y-e i +u) + ^ (60) 

i<N 

This expression has practical advantages: There is no need of truncation in the single-particle 
space, as far as the Green's function is properly calculated. Furthermore, the boundary 
condition imposed on the Green's function provides the exact treatment of the continuum 
states [31, 32]. Normally, the energy of the occupied orbital is negative, e% < [i = 1, • • • , N). 
Thus, the Green's function in the second term in Eq. (60) always has a damped behavior 
in an asymptotic region (r — > oo) because of its negative argument £j — uj < 0. However, the 
first term may have an oscillatory behavior for + oj > 0, which is provided by the outgoing 
boundary condition in . 

An impulsive external potential associated with the function F(f), 

Vl (x) = -sF(f)5{t), (61) 

produces the density response p\{x) as Eq. (47). Here, the parameter s has the dimension 
of ML 2 T^ 1 [F{f f )\^ 1 . The following quantity measures the collectivity of the response: 

R F (t) = ^- J d 3 rF(f) Pl (f,t). (62) 

The Fourier transform of Rp(t) is given by 

R f (oj) = J R F {t)e iut dt = j d 3 r J d 3 r'F(r)U(f,P;uj)F(r / ). (63) 
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Assuming u > and using the relation (x + in) 
is obtained from the imaginary part of Rf(uj). 

S F ( U ) = J2MF\0)\ 2 S(u-(E, 

n 

3.3. Real-time method 
According to the response function formalism in Sec. 3.2, we need to construct the density- 
density response function II by solving the Dyson-type equation (51). The required numerical 
task significantly increases as the dimension of the response function l~I(r, r'; oj) increases, and 
it has been practically prohibited for non-spherical systems. In contrast, the real-time method 
solving the TDKS equation in real time provides a practical and efficient tool for calculation 
of the strength function [33-36] . The method is based on the numerical integration of the time 
evolution of the Kohn-Sham orbitals, described by the TDKS equation (28). The external 
perturbation v(f,t) is taken to be small enough to validate the linear approximation. A 
good account of the continuum is given by the use of complex absorbing potential [34, 35]. 
Recently, the canonical-basis real-time method has been developed for open-shell nuclei with 
the BCS-like pairing. This is based on the diagonal approximation for the time-dependent 
pair potential [25]. In this paper, we cannot present the details of these methods, but readers 
should be referred to Refs. [25, 34, 35]. 

To calculate the density-density response function II, one must evaluate the residual kernel 
w(f,r') = 5v ks (f)/5p(r') (Eq. (53) in the adiabatic approximation), which is a tedious task 
for realistic nuclear EDFs. In real-time method, all we need to evaluate is the KS potential 
Vks[p{t)]- This is a practical advantage in the real-time method. However, the real-time 
method often encounters a problem in numerical stability. This is especially serious in nuclear 
physics, because there is no static external potential to hold the nuclear center of mass at 
a fixed position. Since the translational motion has no restoring force, the moving system 
eventually hits the boundary of the space and produces a spurious contribution to physical 
quantities. Therefore, it is desirable to develop a practical methodology in the frequency 
domain (representation), keeping the advantageous features in the real-time method. This 
is the finite amplitude method (FAM) [39], which is presented in Sec. 3.5. 
To illustrate the basic idea of the FAM, in the followings, we recapitulate the standard 
density matrix formulation and its particle-hole representation. 

3.4- Matrix formulation in the particle-hole representation 

The most standard formulation of the density response is a matrix formulation [19, 27]. We 
start from the TDKS equation (32), where h s (t) contains an external perturbation v\(t). 
Provided that vi(t) is weak, we may linearize the TDKS equation with respect to v\(t) and 
to the density response pi(t). 

p(t) = p + Pi(t), (65) 

h s (t) = h + v sc {(t), (66) 

where ho = h s [po] is the static KS Hamiltonian at the ground state and v sc f(t) is a self- 
consistent effective field induced by density fluctuations, Eq. (55): 

v sc t(t) = ui(f) + v kS:1 (t), (67) 



= Vx — iirS(x), the strength function 
-Eq)) = —ImR F (aj). (64) 
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where t>ks,i(i) = w(t) • pi(t) = 8v]<_ s /8p ■ pi(t). It should be noted that Vks,i(£) has a linear 
dependence on pi(t). Using the stationary condition of the ground-state density, [ho, po] = 0, 
this leads to a time-dependent linear-response equation with an external field, 

^Pi(i) = [ho,Pi(t)\ + [v sc {(t),p ], (68) 
which can be written in the frequency domain as 

u pi(io) = [h ,pi(uj)] + [v sc {(u),p ]. (69) 
Here, we decompose p\{t) and v sc f(t) into those with fixed frequencies: 

Pi(t) = ^{wi^je^ + ^VlHe^}, (70) 

v S cdt) = ^{r/t; scf (a;)e-^ + r / * ? ; s t cf ( W ) e ^}, (71) 

where we have introduced a small dimensionless parameter r\. v sc {(oj) is a sum of v\(oj) and 
^ks,i(w) = Svy^/Sp- pi(u). Note that the transition density pi(w), the external field v\(uj), 
and the induced field t>k s ,i(w), are not necessarily Hermitian in the w-representation. 
The time-dependent KS orbitals as solutions of Eq. (31) are written as \ipi(t)) = + 
\ipi t i(t)), where \(j)i) are time-independent eigenstates of the ground-state KS Hamiltonian 
ho\<Pk) = ffcl^fc)- A proper gauge choice should be adopted to make the stationary eigen- 
states \(f>i) consistent with d\ipi)/dt = 0; e.g., ^ = e^-. Then, the time-dependent density 
response is 

N 

pi(t) = P {t) - po = E ilMWhl + \<k)((Mt)\} ■ (72) 
i=i 

are expanded in the Fourier components as 

\Mt)) = E {«M)e- iw * + V^Y^e^} . (73) 

Ul 

and the density response at the frequency u> is given by 

piM = E + l&X^MI}. ( 74 ) 

The orthonormalization of the TDKS orbitals leads to the fact that only the particle- 
hole (ph) and hole-particle (hp) matrix elements of p\{oo) are non-zero. Namely, (pi)ij = 
(pi)mn = for i,j < N and m,n> N. Thus, without losing generality, we can assume that 
the amplitudes, \Xi(oj)) and \Yi(u)), can be expanded in the particle orbitals only; 

\Xi(u)) = E \<f>m)X mi (u;), \Yi(u)) = E \MY£i(u). (75) 

m>N m>N 

Using the M x N matrix (p a i = (a\(f)i) for the hole orbitals, and the M x (M — N) matrix 
for the particle orbitals (p am = (a\4> m ), the matrix p\{oS) can be expressed by 

pi (w) = tpX(p ] + (pY T ip ] . (76) 

From this expression, it is apparent that the ph and hp matrix elements of p\{oo) represent 
X(uo) and Y(oj), respectively. 
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If we take ph and hp matrix elements of Eq. (69), we can derive the well-known linear 
response equation in the matrix form [19]; 

?{('-)-(i-JU^)-teO- m 

Here, the matrices, A and B, and the vectors, / and g, are defined by 

Ami,nj = (^m — ^i)^mn^ij ~l~ ^mj,ini ^mi,nj = ^mn,iji (78) 

fmi(uj) = (</>m\vi(w)\(f>i), g m i(uj) = (<pi \vi (bj) \ <j) m ) ■ (79) 

The residual kernel w is often called Landau-Migdal interaction, defined by 



Wmk,il = (<Pr, 



dpik 



\<k). (80) 

P=Po 



In nuclear physics, this matrix formulation is also known as the random-phase approximation 

(RPA). The matrix S = in Eq. (77) is Hermitian and called "stability matrix" 

because its eigenvalues characterize the stability of the ground state determined by the 
solution of the KS(B) equations. If S is positive definite, the ground state is stable, thus 
corresponds to a (local) minimum. In contrast, if S has negative eigenvalues, the ground 
state with po is not a minimum, and there exists another true ground state. 
In practical applications, the most tedious part is calculation of the residual kernel, w m j n i 
(w m n,ij)- These elements are two-body-type matrix elements with four indices. Their calcu- 
lation is often the most demanding part in numerical calculations. In the next subsection, we 
propose an alternative approach to a solution of the linear-response equation (69), without 
the explicit evaluation of the residual kernel. 

3.5. Finite amplitude method 

Let us remind ourselves that Eq. (77) was obtained by expanding v sc f(uj) with respect to 
X m i(oj) and Y m i(u). The essential idea of the finite amplitude method (FAM) is to perform 
this expansion implicitly in the numerical calculation. 
Equation (77) reads 

(e m -€i- uj)X mi (uj) + (v ksA ) mi (uj) = -(vi) mi (oj), 

(81) 

(e m - £j + U})Y mi (u) + (t>ks,l)«n( w ) = -(fi)imM, 

where fk s ,i( w ) = d v ks/$P ■ Pi{^)- I n the FAM, instead of performing the explicit expansion 
of fks,i(w) with respect to X and Y, we resort to the numerical linearization. Now, let us 
explain how to achieve it. 

For given amplitudes X and Y, t>k s ,i( w ) can be numerically calculated using the finite 
difference with respect to a small but finite real parameter m 

Vks,i(u) = - (v^lp^uj)] - v^Ipo]) , (82) 

where p v (oj) = po + r}p\(oj). The parameter r\ should be chosen small so that the second and 
higher-order terms in fks[/?r?] are negligible. In the first order in w, p v can be expressed by 
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p v (u) = ip' v ipl, where 

fi, = <p + Ti<pX(u>), 4>l = (<p + W Y*(uj))i (83) 

Namely, fksl/ 9 ??] m Eq. (82) is evaluated simply by replacing the ket states \<f>i) with \<j>i) + 
vYjm \4>m)X mi (uj), and the bra states (<^| with (<^| + r]Y, m ( < Prn\Y m i(^)- Regarding the KS 
potential v ks as the functional of (p and <p', Eq. (82) is rewritten as 

Vks,iH = - (tfks - v ks [p,^}^ . (84) 

The most advantageous feature of the present approach is that it only requires calculation 
of the KS potential, fksfp^]- This should be included in the computer programs of the static 
DFT calculations. Only extra effort necessary is to estimate the KS potential with different 
bra and ket single-particle states, tp' and tp'. Therefore, a minor modification of the static 
DFT computer code will provide a numerical solution of the linear density response. This is 
the essence of the FAM. 

Using these numerical differentiation, both sides of Eq. (81) can be easily obtained by 
calculating the ph and hp matrix elements of the KS potential v ks . Since these are inhomoge- 
neous linear equations with respect to \Xi(uj)) and (Yi(uj)\. we can employ a well-established 
iterative method for their solutions. See Sec. 5.1.3 for more details. 

A typical numerical procedure is as follows: (1) Fix the frequency ui and assume initial 
vectors (n = 0), X^(u) and y£)(u>). (2) Update the vectors, X^ +1) (u;) and Y^ +1 \oj), 
using the algorithm of an iterative method. (3) Calculate the residual of Eq. (81). If its 
magnitude is smaller than a given accuracy, stop the iteration. Otherwise, go back to the 
step (2). 

To calculate the strength function with respect to the Hermitian operator F, we should adopt 
vi(oj) = F. After reaching the solution (X m i,Y m i), the strength function (64) is obtained by 

S F (u) = —ImR F (u), (85) 

7T 

where 

R F (u>) = tr[Fpi(w)] = £ {F m X m (io) + F mi Y mi (uj)} . (86) 

mi 

3.6. Quasiparticle formalism with pairing correlations and FAM 
In case that the pairing correlations play essential roles, we extend the previous TDKS 
formalism to the TDKSB formalism in Sec. 2.7. This is straightforward. Starting from the 
equation for the time-dependent generalized density matrix, Eq. (44), we follow the same 
procedure as that in Sec. 3.4. In this section, all the quantities must be defined in the body- 
fixed frame rotating in the gauge space, which are expressed with the prime (') in Sec. 2.7 
6 . We omit the primes here for simplicity. 



This corresponds to the "moving-frame harmonic equation" in Sec. 4.5 
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The TDKSB Hamiltonian % s contains the unperturbed % s , a time-dependent external 
perturbation Vi(t), and induced field H s ,i(t). 



Hs[R](t)=H° 9 +H a ,i(t) + V 1 (t). 



(87) 



This leads to the generalized density matrix R(t) = Rq + R\(t), where Ro is the ground-state 
density in the rotating frame. Following the same arguments as in Sec. 3.4, we may derive 
the linearized TDKSB equation for the generalized density response R±(lo), 



uR±{u)= [HPs-^Riiu)] +[H s ,i{u) + V 1 (uj),R ]. 
Using the matrix notation of Eq. (36), the qp orbitals are expressed as 7 

= * + ¥(t) = $ + *i(t), 



(88) 



(89) 



with R = = 1 - <S>&. 

can be expanded only in terms of the "unoccupied" 
("occupied") static orbitals <& ($), thus written as 



^(t) = $Y1 { X (u)e- iwt + Y*(oj)e iujt } 
^(t) = {X*{u)e +iljt + Y^e-^} . 



(90) 



X vll t{u) and Y uu i(oj) are M x M matrices, which must be anti-symmetric because of the 
unitarity of W(t) = The density R(t) = V(t)&(t) is expanded up to the first 

order in ^i, which gives Ri(t) = *S>i(t)& + $ty\(t). Substituting Eq. (90) into this, we have 



Ri(u) = + $y t &. 



(91) 



From this expression, one can see that the "unoccupied" -"occupied" matrix elements of 
R\{oj) are expressed by X, and the "occupied" -"unoccupied" matrix elements are given 
by Y. This is analogous to Eq. (76). Using the unitarity of W = (<&,<£) and the following 
relations, 

(H° s - MA) = (H° s - fjJ\T) § v = -E u $ u , 

then, Eq. (88) leads to the linear density response equations: 

(E v + E v , - u)X vu , + (H S ,!)ZH = -(Vi)^(w), 
(E v + E v , + U )Y W , + {U s ,i)fA") = -(ViC(w), 



(92) 



(93) 



where 



\20 



- vv' 



(Vi 



\20 



$ f Vi(w)$ 



vv' 



(94) 
(95) 



If we expand % 2 ®i and in terms of X and 1", we reach the familiar expression of the 
matrix form, similar to Eq. (77). The A and B matrices are given by the qp energy and 



7 Here, we assume a proper choice for the gauge parameter S„„/ to make a stationary solution <3> 
time-independent. 
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the residual kernels, 

where the "unoccupied" -"occupied" and "occupied" -"unoccupied" components of the gen- 
eralized density, it! 20 and i? 02 , are defined by 



o20 



R 02 



0,02 



(97) 

The finite amplitude method (FAM) for the qp density response is presented in Ref. [40]. 
Here, we recapitulate the essential idea and the result. Instead of calculating the residual 
kernels in Eq. (96), (% S) i) 20 and (% Si i) 02 in Eq. (93) are numerically obtained by the finite 
difference. First, we define the 77-density R n {oo) as 

R^oj) = R + vRiico) = %(cu)^l(u), (98) 

where 

^( w ) = ® + v <f>X(oj), ft(w) = (fc + f/Sr^w))* . (99) 
Then, the induced residual fields are given by the following FAM formula: 

^20 = ^ Ks[R v ]-Ks[Ro} ^^ h w = _Q\ H a [Rr,]-H a [Ro] ^ (10Q) 

Equivalently, the FAM formula can be written in terms of the qp orbitals as 

«5 = «Wti -«•"■ *'U (101) 
tt w.[$;,$t]-w.i*,*\ (102) 

77 

A computer program for stationary solutions of the KSB equation is able to construct the 
KSB Hamiltonian "H s [i?] from the qp orbitals <I>. Thus, a small extension of the code to 
calculate % s for different <!> and $t allows us to turn the static KSB code into the one 
for the linear response calculation. The FAM significantly reduces the programming task of 
developing a new code [40, 41]. It turns out to save the enormous computational resources 
as well, in linear-response calculations for deformed nuclei [41]. 

4. Beyond perturbative regime: Large amplitude dynamics 

Nuclei exhibit a variety of collective phenomena with large-amplitude and anharmonic nature 
in the low-energy region. For instance, the nuclear fission is a typical example for such a large- 
amplitude collective motion, that a single nucleus is split into two or more smaller nuclei. To 
describe these large-amplitude phenomena, we are aiming at developing a practical theory 
to extract a few optimal canonical variables, to describe the slow collective motion, which 
are well decoupled from the other fast intrinsic degrees of freedom. Then, upon the obtained 
submanifold, the collective Hamiltonian is constructed with microscopic determination of 
the collective mass parameters and potentials, to calculate observables in nuclear collective 
phenomena. 

There have been extensive efforts in nuclear theory for such purposes (See recent review 
papers [42, 43]). In this article, we present a classical theory of the large amplitude collective 
motion. The contents in this section is mostly based on former works [42, 44-46]. 
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4-1. Basic concepts of a decoupled collective submanifold 

As is shown in Appendix of Ref. [42], the TDKS(B) equations are identical to the classical 
Hamilton's equations of motion with the canonical variables 7T a }a=i,- ,N C - The number of 
independent variables iV c are in the order of M 2 for the description of the TDKSB dynamics. 
Since M is in principle infinite without the truncation, N c could be huge for description of the 
large-amplitude motion. Thus, it is desirable to extract a few canonical variables which are 
approximately decoupled from the other degrees of freedom. These variables are supposed to 
describe decoupled collective motion of the many-body system. There are several equivalent 
ways to present the basic concepts and formulation of the theory. 

We assume the collective motion of interest is a slow motion which allows us to truncate 
the classical Hamiltonian under the expansion with respect to momenta. Up to the second 
order in momenta 7r, the system is described by the Hamiltonian 

H(£,, r ) = ±B a f ) (Z)wf l + V(Z). (103) 

The summation with respect to the repeated symbol for upper and lower indices is assumed, 
hereafter. The reciprocal mass tensor B al3 is defined by 

Ba0 9 2 ff(g,7T) 

The mass tensor B a p is defined by B a pB^ = (52, as the inverse matrix of B a/3 . We are trying 
to find a collective submanifold present in the classical Hamilton system described by H in 
the form of Eq. (103). 

4-1-1. Point transformation. In general, The main aim of the theory is to find the 
canonical transformation 

teW-K^M. (105) 

where the {q 11 } are assumed to be divided into two subsets, q l , i = 1, ■ ■ ■ , K and the rest q a , 
a = K + 1, • • • , N c , which are decoupled with each other. Namely, if the system is located at 
q a = and q a = at time t = 0, then the time evolution should keep q a (t) = 0. Of course, in 
reality, the decoupling is not exact. We want to find an approximately decoupled manifold. 
First, we limit ourselves to the point transformations. 

?" = /"(£), e = 9 a {q)- (106) 

In the point transformation, the transformation for conjugate momenta are given by 
derivatives of the functions / M and g a . 

Pn = 9 a ^a, TTa = ffcPp, (107) 

where the comma indicates a partial derivative, g a ^ = dg a /dq^. The canonicity is guaranteed 
by the conservation of the Poisson brackets, which is easily proven by using the chain-rule 



(104) 
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relations: 

= f%9*v = %- (108) 

Substituting the point transformation of Eq. (106) into Eq. (103), the Hamiltonian in the 
new variables becomes 

H(q,p) = X -B^{q) VixVv + V(q). (109) 
The reciprocal mass parameter transforms like a tensor of the second rank. 

B^(q) = ]%B^jy (110) 

4-1-2. Decoupling condition under a point transformation. The decoupling condition is 
that, if the system is located on the collective submanifold (q a = p a = 0), it stays within the 
submanifold, namely, q a = p a = 0. From the Hamilton's equations of motion derived from 
Eq. (109), 

r) fT 

■a fjai^ I Dab„ 

P - (HI) 
■ a _ _ _y _ i Dja' 

dq a ~ ,a 2 ' a M 
we have the following conditions for the decoupling: 

B ai = 0, (112) 
V a = 0, (113) 
B% = 0. (114) 

The first condition, Eq. (112), tells us that the reciprocal mass tensor must be block diagonal 
and has no coupling between the collective space (q 1 with i = 1,- ■ ■ , K) and the intrinsic 
space (q a with a = K + 1, • • • , N c ). The remaining two conditions comes from the absence 
of the force perpendicular to the collective surface. The conditions for the mass tensor, Eqs. 
(112) and (114), imply that the decoupled submanifold is geodesic with the metric given by 
the mass tensor By. Namely, the following quantity 

Ax = J \f~Bdq 1 A • • • A dq K , B = det(Bij), (115) 

with a fixed boundary 8Sk has a minimum value, SAk = 0. See Ref. [42] for the proof. 
Utilizing the chain rule, the force condition, Eq. (113), can be rewritten as 

V a = V tli f% = %f] a . (116) 

This is the condition obtained in the zero-th order in momenta. In the one-dimensional case 
[K = 1), it is 

V a - \if^ a = 0, (117) 

where Ai = dV/dq 1 . This is nothing but the minimization of the potential V(£,) with a 
constraint on the collective coordinate q 1 = f 1 . 

s{v(0-\ifHO}=o- (us) 
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4-2. Local harmonic equations (LEE) 

If all of the decoupling conditions, Eqs. (112), (113), and (114), are satisfied, it provides an 
exactly decoupled collective submanifold. However, in realistic situation, the exact decou- 
pling is not realized except for trivial collective motions, such as the translational motion. We 
are more interested in situations with approximate decoupling. Among the three decoupling 
conditions, Eq. (114) comes from the coefficients in the second order in momenta. Here, we 
build the theory by ignoring this second-order condition, based on the mass condition (112) 
and the force condition (113). 

Let us start from the chain-rule about the derivative of the potential, 

V a = V u f% = V a g%, (119) 

This indicates that the first derivative has a property of the covariant vector. However, the 
second derivatives are known to be not a tensor with respect to the general point transforma- 
tion. As is well-known in the general relativity, we should introduce the covariant derivative, 
to keep the tensorial property. The covariant derivative is defined by 

V ** 55 dV ' a dJ V ' a = V <# ~ (120) 
using the parallel transport of the vector V a (£) for £ — > £ + d£, 

sv a = r^vpdp. (i2i) 



In order to make the covariant derivatives V a p a tensor of the second rank, the affine 
connection must follow the transformation: 

n p = f%g^gy% + f%9% P i ( 122 ) 

Now, we assume that the coordinate system {q^} is geodesic, namely, r„ p = 0. This leads 
to the affine connection 

and the covariant derivatives 

V a p ee V a p - Tl p V, = V a p - r a pV„ % v = V, u , (125) 
which can be even simplified because of Eq. (113), as 

V a p = V aP - r tap Vj. (126) 
Since these covariant derivatives are tensor, they must transform as 

v a( 3 = WM- ( 127 ) 

Multiplying the reciprocal mass tensor, we have 

M% = B^V nP = Bivv^&tt = Mtglf% (128) 

from which we easily obtain the following equations. 

My% = M$f% M%^ = Mlg% (129) 

Now, let us use the decoupling conditions, Eqs. (112) and (113). Taking u = i (collective 

coordinate) in Eq. (129), the decoupling conditions tell us that the matrix My = B^ p V pu 
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is also block diagonal, M. l a = and M.® = 0. Therefore, we reach the following equations, 
which we call "local harmonic equations" (LHE). 

MV; n M'.r r M%fi = M\g a j. (130) 

In the case of K = 1, it is written as 

M%f\ = u 2 f% M%g? = u 2 g% (131) 

where the frequency is given by uj 2 = ftA\. The solution of the LHE provides a tangent vector 
of the collective submanifold, f ' l a and g°-. 

The LHE generalizes the secular equation of the harmonic approximation around the 
potential minimum to that at non-equilibrium points. At the equilibrium (V a = V jfl = 0), 
it automatically reduces to the normal harmonic approximation, because the covariant 
derivatives become identical to the second derivatives at the equilibrium, V a p = V a p. 

4-2.1. Practical solution of LHE. To solve the LHE (130), we need to calculate the affine 
connection of Eq. (124), which contains the curvature f l aj3 , in the covariant derivative V a p. 
Since the solution of the LHE provides only f l a and this cannot be given by the LHE 
itself. However, the curvature f l a g can be eliminated in the following procedure [45]. Here, 
let us discuss the case of K = 1, for simplicity. In this case, we have a single collective 
coordinate q 1 . We take the derivative of B 11 = B a/3 f) a f \ with respect to q 1 . 

B]l = Bj/i/M + 2B«y]j^ v (132) 

Using this equation, we may eliminate the curvature terms in the LHE (131). The LHE can 
be rewritten in the same form as Eq. (131), but Mp and uj 2 can be replaced by 

M% ee Bnv aP + \b^v, 

1 (133) 
J = B ^v u + -BfV^. 

In this way, we can eliminate the curvature terms. The price to pay is the calculation of 
B°p . The eigenfrequency uj 2 is obtained by solving the eigenvalue equation (131). Thus, We 
do not need to calculate B 1 } . 

4-2.2. Riemannian connection. In this article, we adopt the affine connection of Eq. (124), 
which assumes that the decoupled coordinates is geodesic (T„ p = 0). Instead of 

Eq. (124), the Riemannian connection may be adopted, in a similar manner to the general 
relativity. The Riemannian connection is given in terms of the metric tensor B a p as 

T % = 2 BaS ( B /3S,i + B -yS,p ~ #87,5) • (134) 

In Ref. [42], this was discussed in details. In fact, the decoupling conditions, Eqs. (112), (113), 
and (114), may lead to the LHE identical to Eq. (130) with the connection Pg 7 replaced 
by Eq. (134). However, the Riemannian formulation of the LHE has a problem for the case 
that the Nambu-Goldstone modes exist [44], which will be discussed in Sec. 4.3 Thus, in the 
followings, we focus our discussion on the LHE with the affine connection of Eq. (124). 
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4-3. Treatment of constants of motion: Nambu-Goldstone modes 

Since nuclei are self-bound system without an external potential, the nuclear DFT provides 
a ground-state density distribution which spontaneously violates the symmetry, such as 
translational and rotational symmetries. The spontaneous breaking of symmetry produces 
the Nambu-Goldstone (NG) modes which correspond to trivial (spurious) collective degrees 
of freedom. Therefore, we are mostly interested in the extraction of the collective degrees 
of freedom which are separated from (perpendicular to) these NG modes. In this section, 
we show that the LHE presented in Sec. 4.2 properly separate the NG degrees of freedom 
from other degrees of freedom. However, to achieve this, we need to lift the restriction to the 
point transformation and extend the point transformation to allow the second-order terms 
in momenta [44]. 

4-3.1. Extended adiabatic transformation. The restriction to the point transformation is 
lifted by expansion with respect to momenta it. Equations (106) are generalized by 

t = P{t) + \f (l) ^{t> a itp + 0(^), (135) 

r = 9 a (q) + \g {l) ^ u {q)P»Pv + 0{p A ). (136) 

The transformation of the momenta is given by Eq. (107), since the terms cubic in momenta 
do not play a role in the modification of the theory. Using Eq. (107), the independence of 
the variables, d^ a /dtrp = 0, requires the relation 

gW^g^ = -fW&gfr. (1 37 ) 

From the canonicity condition q u }pB = 0, we also find 

= ^W, ( 138) 

The Hamiltonian (103) is transformed to, up to second order in p, 

H(q,p) = V(q) + \b^ v , B^ = f%B^f^ + Vj 1 ^. (139) 

The major difference between the use of the extended adiabatic transformation and a point 
transformation is the modification of mass parameter, 

B afi = g^B^g^ = B a/3 -V^f {l) ^ (140) 
= B a P-Vif mia P. (141) 

Here we have used Eqs. (137) and (113). The LHE has the same form as the Eq. (130), after 
replacing B al3 by B a ^ . 

4-3.2. Constants of motion; cyclic variables. Suppose a classical variable V(^,tt), which 
correspond to one-body Hermitian operators P in the quantum mechanics, is a constant of 
motion. In the followings, the conserved quantities are classified into two categories. Adopting 
the classical canonical variables in Ref. [44], if P has real matrix elements in the qp basis, 
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V can be expanded as 

V(M =P {0) (0 + ^ (2)Q V^ + O(^ 4 ). (142) 

On the other hand, if P has imaginary matrix elements, 

V(t,ir)=V^ a K a + 0(^). (143) 

The conservation of V indicates that the Poisson bracket between V and H should vanish. 
From this, terms of the zeroth and first order in tt give the following identities. 

V^ a V a = 0, (144) 
V (0) B al3 _ <p(2)aPy a = q ( 145 ) 

The equations (144) and (145) hold at arbitrary points in the configuration space. 
We assume that the variables describing these constants of motion correspond to the canoni- 
cal variables (q 1 ,pi). The collective variables (q l ,Pi) of interest are supposed to be orthogonal 
to both these variables and the intrinsic variables (q a ,p a )- Thus, we divide the set {q^,Pfj,} 
(fx = 1, ■ ■ ■ , N c ) into three subsets, the collective coordinates {q l ,Pi}, i = 1, • • • , K, the cyclic 
coordinates {q^pi}, I = K + I, ■ ■ ■ ,K + M, and the non-collective coordinates {q a ,p a }, 
a = K + M + ,N C . 

In nuclear physics applications, we are often interested in the large amplitude collective 
motion at a given value of q 1 , such as the given total angular momentum, q 1 = (J x ) = J, 
and the given number of particles in the superfluid systems, q 1 = (N) = N. In this case, Eq. 
(116) should be modified with additional constraints with respect to q 1 as 

V, a -Vtf a -Vrf a = Q. (146) 

Here, f J a are given by the 2qp matrix elements of the symmetry operator. The non-trivial 
collective coordinates of interest, f l a , are determined by the solution of the LHE (130) with 
the reciprocal mass tensor of Eq. (141). We should solve Eqs. (146) and (130) self-consistently. 

4-3.3. Separation of cyclic variables as zero modes. Now, let us prove that ff a (g a j) pro- 
vides the zero-frequency solution (ui = 0) for the LHE. We start from the case that the 
coordinates q 1 are conserved, with / J (£) =p(°)(£) and fW Ia ^) = V 1 - 2 ^ {£) . This cor- 
responds to the case of most practical interests in nuclear physics, such as the angular 
momentum and particle number. 

M a f i a = B^V^ a = {B a y% - V^fW^fijV^ 

= (B a y* a - v^ Ia y^v. lP = {B a yi a - v a f^)v^ = o. 1 

Here, Eq. (138) was used in the third equation, and Eq. (145) was in the last equation. Thus, 
automatically becomes a solution of the LHE with oj = 0. 
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Next, we discuss the case that the momenta pi are conserved, with <?^(£) =V^ a {^). 
Differentiating the chain relation <7^/^ = 6p with respect to q u , we obtain 

9%r, = ~(M,f%- (148) 

Differentiating Eq. (144), we have 

VapgJ + V, a g%f% = 0. (149) 
Utilizing these equation, we may prove 

~ ~ (150) 

= B^iV^ - f^Vsg 5 ^) = B^iV^g? + f»g%V s ) = 0. 

Therefore, the gf are zero-frequency solutions of the LHE. 

The separation of the NG modes are guaranteed in the LHE with the covariant derivatives 
V a p of Eq. (126) and the reciprocal mass tensor B a/3 of Eq. (141). 

4-4- Gauge invariance 

The basic formulation to determine the collective submanifold is given by Eqs. (116) and 
(130). In the case of the one-dimensional collective coordinate (K = 1), these equations 
provide a unique solution, except for the scale of the collective coordinate q 1 . However, for 
the multi-dimensional collective manifold (K > 0), the solution of Eqs. (116) and (130) are 
not unique. In fact, there is a gauge invariance similar to what we observed in Eqs. (31) 
and (38). For a pair of collective variables (q k ,Pk) and (q l ,Pi), k ^ I, we may adopt a point 
transformation 

q k ->q k + cq l , p, -»• p x - cp k , (151) 

with an arbitrary gauge parameter c, keeping the other variables unchanged. Let us show the 
transformation of Eq. (151) keeps the formulation of Eqs. (116) and (130) invariant. Namely, 
(q k + cq l ,pk) and (q l ,Pi — cpk) instead of (q k ,Pk) and (q l ,pi) also provides a self-consistent 
solution for Eqs. (116) and (130). 

Since the transformation (151) gives g? — > g°i — cg a k , the derivative of the potential, Vj = 
V tQ .g", is transformed as Vj — > Vj — cV^k- From this, we can immediately see the invari- 
ance of Eq. (116), using the transformation f k a — > f k a + c/' a . The matrix Ai 1 ^ in the left-hand 
side of Eq. (130) is also invariant under Eq. (151). In fact, both B a/3 and V a p are separately 
invariant. In contrast, the matrix M 1 - in the right-hand side of Eq. (130) transforms as 

M k -> M k + cM\, Mj^ Mj - cM{ , M{ -»• M\ , 
Mf -> Mf + cM\ - cM\ - c 2 M{ 

for i ^ / and j / k. This can be easily obtained from the relation M 1 - = f^M^g^. These 
relations prove that {P a }i=i,~ ,k with f k a replaced by f k a + cf l a also provides a solution of 
Eq. (130). In the same manner, we can prove that {<7"}i=i, -- ,k with gj replaced by g°\ — cg a k 
is a solution as well. This gauge invariance is present for any pair of collective variables (k, I), 
thus for an arbitrary linear point transformation. 

In the case that the cyclic variables (q 1 ,pi) exist, the gauge invariance is present even for 
K = 1. Suppose q 1 is a collective coordinate, which is a self-consistent solution of Eqs. (146) 
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and (130) with the mass tensor of B a/3 . Then, the following transformation provides another 
solution: 

q 1 — > q 1 + cq 1 , pj ->■ pi — cp\. (153) 
The proof is given by exactly the same argument done for Eq. (151). 

This gauge invariant property tells us that we need to fix the gauge parameter c. For instance, 
a possible choice could be requiring Vu = which was adopted in Ref. [46]. One can make 
other choices if they are more convenient [46-48], and the physical quantities should not 
depend on this choice. 

4-5. Moving-frame harmonic equation (MFHE) 

Let us summarize the formulation we obtained so far. The present formulation can be 
regarded as the harmonic equations with the moving-frame Hamiltonian 

H M {£, tt) = H{£, tt) - Xiq 1 - \ iq \ (154) 

Equations (146) and (130) can be rewritten as 

S{H M } n=0 = 0, (155) 

(M M )P: a = (M M ) i j ff ) , or (M M r^ = {M M ))g% (156) 

Here, the matrix (Mm)'q is a product of the mass and potential, given in the same way as 
Eq. (103) but with H M - 

B M= B ^ = -^Jt' V M (£) = H M {£,* = Q). (157) 

It turns out that the LHE becomes identical to the harmonic equation at the equilibrium with 
Hm- Therefore, we may call this formulation "moving-frame harmonic equation" (MFHE). 
It should be noted that the terms —Xiq 1 — Xiq 1 are not merely the constraints. These terms 
changes the mass parameters and the potential. The theory of the MFHE is basically equiv- 
alent to the gauge-invariant formulation of the adiabatic self-consistent collective coordinate 
(ASCC) method [46]. 

From this moving-frame formulation, it is evident why we use the Hamiltonian in the rotat- 
ing frame, H — XN, in Sec. 3.6. The same argument is also applicable to the quasiparticle 
random-phase approximation (QRPA) in the superfluid phase [19, 27]. Since the ground 
state does not correspond to the equilibrium of the energy surface (dE/dN / 0), the QRPA 
is a harmonic approximation at a non-equilibrium state. According to the present theory, 
requirements of the covariance and the extension of the point transformation defines the 
moving frame in which the QRPA should be formulated. 

4-5.1. Practical solution of MFHE. The theory to define a decoupled submanifold con- 
sists of Eqs. (155) and (156): The first equation (155) is the potential minimization with 
constraints on q l and q 1 , which defines the position £. The second equation (156) defines the 
normal modes p a (g°) at the same position £, which should provide p a used in Eq. (155). 
Therefore, these equations should be self-consistently solved. 

Let us discuss the K = 1 case in more details, how to construct the MFHE matrix (Mm)^- 
The MFHE (156) contains higher-order terms which are not present in the LHE discussed 
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in Sec. 4.2.1: fj afi , ft 1 ) 1 "?, f a/3 , and fW la P. Among these quantities, f afi and f^ Ia/3 are 
calculable if we know the operators corresponding to q 1 explicitly, such as the particle number 
and the angular momentum. The curvature f 1 ^ can be eliminated by the same procedure 
as that in Sec. 4.2.1. Thus, the remaining unknown quantity is /( 1 ) lct ' 3 . 
Although we do not have a general principle to determine fW la ^ l there may be possible 
prescriptions. In the case that there is a single constant of motion q 1 , the canonicity condition 
of Eq. (138) gives constraints whose number is same as the number of the index a, namely 
the number of 2qp states. Using these constraints, possible prescriptions are, for instance, 

(1) Diagonal assumption: Assuming fi 1 ) 1 ^ = f( l ) la ^ a P j fW la can kg determined by 
Eq. (138). 

(2) Strong canonicity condition: Both g 1 and q 1 are assumed to be represented by one- 
body operators Q l and Q 1 , respectively, where Q 1 is explicitly known. Then, requesting 
[Q l ,Q r ] = can determine fi 1 ) 1 ^. 

In the numerical applications in Sec. 6, we adopt the prescription 2 to examine the effect of 

Effect of this term turns out to be negligibly small for the multi-0(4) model [46]. 
After eliminating the curvature terms, the MFHE can be rewritten in the same form as Eq. 
(131), with 

(M M )$ = B^(V nP - Vjf^) + \&gV a 

(158) 

(u M f = B ll V u + -SyVx. 

The equations equivalent to these have been solved in Refs. [46-48], with the second 
prescription given above. 

5. Giant resonances studied with Skyrme EDFs in the linear regime 

Applications of the TDDFT have been mostly studied in the linear response regime. In this 
section, we show selected results of the applications of the Green's function method (Sec. 3.2) 
and the finite amplitude method (Sec. 3.5) for nuclei without the pairing correlations, and 
the standard diagonalization method [19] for superfluid nuclei. 

5.1. Giant resonances in the normal phase 
5.1.1. Coordinate- space representation. 

For the Skyrme functionals, which is a functional of local one-body densities, the coordinate- 
space representation is one of the convenient choices [49]. In the followings, we assume f 
involves the spin and isospin indices, if necessary. We adopt the three-dimensional (3D) 
Cartesian grid-space representation in Sec. 5.1.2 adn Sec. 5.1.3. Each KS orbital (pi(r) is 
represented at discretized grid points {xd, yd, Zd). In the linear regime, behaviors of the TDKS 
orbitals Y>j(r, t) in the region far outside of the nucleus are irrelevant in the calculations. This 
is because the density response pi(f,t) vanishes where the KS orbitals in the ground-state 
Mr) = 0: 

= E {\Mm 2 - \Ur)\ 2 } = Efc(r,t)^(r1 + c.c.} = 0. (159) 

i i 
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Thus, we use the 3D grid representation with an 
adaptive mesh [35, 50], to reduce the number of 
grid points in the outer region. See Fig. 2 for such 
an example. 

The forward and backward amplitudes, X m i(cu) 
and Y m i(uj), in the linear response equations (81) 
also possess two indices, (mi). For these, it is con- 
venient to adopt the mixed representation: the 
particle index m > N is replaced by the coordi- 
nate r, but the hole index i < N is kept. The 
number of hole (occupied) orbitals in finite nuclei 
is of order of 100, at most. This mixed represen- 
tation is adopted in the application of the FAM 
in Sec. 5.1.3. 

5.1.2. Application of Green's function method. 
We first show applications of the Green's func- 
tion method. In the case that the KS orbitals are defined in a potential with the spherical 
symmetry v s (r) = vo(r), this is known under the name of "continuum RPA" in the nuclear 
physics [31] and has been extensively utilized to study giant resonances in nuclei [51-53]. It 
should be also noted that the extension to the linear density response in superfluid systems 
has been achieved with the use of the anomalous Green's function [54]. 

For deformed systems, the construction of the Green's function with a proper boundary 
condition involves a significant task [34, 55] and the applications to nuclear systems are still 
very limited. We adopt an approach using a double iterative algorithm [34, 35]. Roughly 
speaking, this is based on the fact that Eqs. (56) and (57) are rewritten in a form of the linear 
algebraic equation with respect to pi(w), In addition, the action of the Green's function, 
\ip(^} = (E) | <p) for a given state \<p) is also given by a solution of the linear equation. 
We solve these linear algebraic equations by using the iterative methods. See Refs. [34] for 
details. 

In Fig. 3, we demonstrate an example of the results of the present iterative algorithms 
for deformed systems. The isoscalar monopole and quadrupole strength functions in 20 Ne 
are calculated with the BKN energy functional that is a simplified version of the Skyrme 
functional [56]. In nuclear binding energy, there is a strong cancellation between the positive 
kinetic energy and the negative potential energy. The large nucleonic kinetic energy plays 
an important role in many phenomena in nuclei. The giant quadrupole resonance (GQR) is 
such an example. Namely, the restoring force for the vibrational motion mainly comes from 
the distortion of the Fermi sphere in the momentum space [57]. 

The GQR shows three peaks in order of K = 0, 1, and 2 in increasing energy (Fig. 3 (b)). 
This is because the ground state has a superdeformed prolate shape with f3 ~ 0.6. The result 
also indicates no low-energy quadrupole vibration except for the NG mode with K = 1. This 
is a characteristic feature of the superdeformation [58, 59]. 

The monopole strength consists of two components: a peak at 15 MeV and a broad hump in 
the energy region of E > 20 MeV. The dotted line indicates the strength of the independent 
particles obtained by Il s . The residual kernel w(u) shifts the two components to opposite 
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Fig. 2 Adaptive grid in the (x, y)- 
plane, used in calculations in the 
following sections. 
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Fig. 3 Strength functions for 20 Ne calculated for complex frequencies ui + with the 
smoothing parameter 7 = 1 MeV [35]. (a) isoscalar monopole strengths. The solid line cor- 
responds to the full response ("TDDFT"), while the dotted line is obtained by neglecting 
the residual kernel, that is indicated by "IPA" (independent particle approximation), (b) 
isoscalar quadrupole strengths. The K = 0, K = 1, and K = 2 quadrupole strengths are 
shown by solid, dash-dotted, and dashed lines, respectively. 

directions. The peak at ui ~ 20 MeV is shifted to lower energy by about 5 MeV. This lowering 
in energy is due to strong coupling to the quadrupole resonance. In fact, the peak lies at 
exactly the same energy as the K = quadrupole resonance (Fig. 3 (b)). 
The calculated single-particle energy of the last occupied orbital is —10.8 MeV. Thus, all 
the high-energy peaks in Fig. 3 are embedded in the continuum. The broad structure of the 
monopole strength function at E > 20 MeV indicates that there is no prominent monopole 
resonance in this nucleus, except for the peak due to the coupling to the GQR. 

5.1.3. Application of FAM. The FAM is a feasible approach to the linear response cal- 
culations with realistic EDF. With a Skyrme-type EDF, the FAM formula (84) tells us to 
calculate the operation of v^ Si i(f,u) in the coordinate space as 

VkB,i(r, w)&(fO = i [h[^*] (r)Mr) ~ , (160) 

with ^Jjj(r) = 4>*(f) + rjY*{u,r) and = 4>i(r) + r}Xi(r,uj). Exchanging the forward 

and backward amplitudes in V^,i(r) and j(r), we may calculate v\ s 1 (r, uj)(pi(f) in the same 
way. Adopting the local external field vi(f,u) = F(f), the strength function is calculated 
from the obtained forward and backward amplitudes, as in Eqs. (85) and (86), 

S F (u) = ~l lm Yl J dr{cp*{r)F\r)Xi{r,uj) + V/(r, w)Ft(r)^(r)} . (161) 

i 

The FAM makes a coding of the linear response calculation much easier than the other 
methods. The FAM does not require explicit construction of the matrix, thus, it significantly 
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Fig. 4 Convergence property of different iterative methods: Relative residue, r n , is shown 
as a function of iteration number n, at complex frequencies of oj = + 0.5i MeV (a) and 
10 + 0.5i MeV (b). See text for details. The figure is taken from Ref. [50]. 

reduces a memory resource requirement. These are the main advantages of the FAM. In 
addition, the computational task scales linearly both with the size of the model space and 
with the particle number. This linear dependence was confirmed in the actual calculations 
as well. Therefore, the FAM may demonstrate its merit for larger systems. 
A disadvantage is the fact that the iterative procedure is difficult to parallelize. Since the 
calculations with different oj are independent, this provides a trivial parallelization with 
respect to oj. This leads to a use of PC cluster systems with 128-256 processors in parallel, 

Choice of iterative algorithms 

To solve the linear response equations (81), an iterative method is utilized. Here, we denote 
this equation symbolically as As = b. For the Skyrme energy functional, the matrix A 
in the r-space grid representation is sparse. Therefore, the iterative methods, such as the 
conjugate gradient (CG) method [60], should work efficiently. However, since we calculate for 
the complex frequency oj, the matrix A is not Hermitian. Therefore, we should adopt one of 
a number of variants of the CG method extended for non-hermitian problems. In Fig. 4, we 
show performance of some of different iterative algorithms: Bi-conjugate gradient (Bi-CG) 
method [60], generalized conjugate residual (GCR) method [61], generalized product-type 
bi-conjugate gradient (GPBi-CG) method [62], Bi-CGSTAB method [63], and Bi-CGSTAB2 
method [64]. The magnitude of the relative residue, 

r n = \b-Ax- n \/\b\ (162) 

is plotted against the number of iterations, for the case of the electric dipole response in 
16 0, The initial vector is taken as x = 0. 

It turns out the convergence property depends on the frequency. At low frequency (oj = 
+ 0.5i MeV), all the solvers except for the Bi-CG method quickly reach the convergence. 
On the other hand, at higher frequency (oj = 10 + 0.5i MeV), only the GCR and the GPBi- 
CG methods lead to the convergence. In Fig. 4, the GCR shows the most stable behavior 
for the convergence, though it requires larger computer memory resources than other meth- 
ods. Recently, we have also tested the generalized product-type bi-conjugate gradient method 
with associated residual (GPBiCG-AR) [65], which indicates a better performance. It should 
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Fig. 5 Calculated (a-c) [50] and experimental (d-f) [67, 68] photoabsorption cross sections in le O, 40 Ca, 
and 2 Mg. We use the SkM* parameter set and 7=1 MeV. 



be noted that we need much smaller number of iteration to reach the convergence in the 
harmonic-oscillator-basis representation [41]. The coordinate space of a relatively large 3D 
box size contains a large number of irrelevant mesh points, which perhaps makes the con- 
vergence very slow. It should be noted that an iterative algorithm based on the Arnoldi 
diagonalization method was proposed for similar problems [66]. 

Nuclear photoabsorption cross sections 

Adopting the electric dipole (El) operator with the El recoil charges as the operator F, we 
calculate the El strength function Sei{E) that is converted into the photoabsorption cross 
section cr a b s (-E) in the dipole approximation. Calculated photoabsorption cross sections for 
spherical nuclei, 16 0, 40 Ca and deformed nucleus 24 Mg, are compared with experimental data 
in Fig. 5. Here, the complex frequency u + 27/2 with the width 7 = 1 MeV is adopted. In 
each nucleus, there is a broad peak in a a ^ s (E) around E = 20 — 25 MeV, which corresponds 
to the giant dipole resonance (GDR). The overall profile of the experimental cross section 
is well reproduced, though the calculated energies of the GDR peaks are underestimated by 
a few MeV. The discrepancy is more prominent for lighter nuclei, which is observed with 
almost all the Skyrme energy functionals [35, 50]. 

For spherical nuclei, the GDR widths calculated with 7 = 1 MeV are narrower than the corre- 
sponding experimental data. This seems to suggest that the spreading width F^, which takes 
account of effects decaying into compound states, such as two-particle-two-hole excitations, 
is slightly larger than 7 = 1 MeV. For the deformed nucleus 24 Mg, the GDR peak splitting 
caused by the ground-state deformation well agree with the experiments, although the mag- 
nitude of the deformation splitting is slightly too large in the calculation. We may interpret 
that the experimental GDR peak around E = 20 MeV is associated with the K = mode, 
and those at E = 22 ~ 25 MeV correspond to the K = 1 mode. The double-peak structure 
of the K = 1 GDR peak is well reproduced. Approximately, the calculated cross section is 
shifted to lower energy from the experimental ones, by about 3 MeV. 
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Fig. 6 Calculated photoabsorption cross sections for (a) 90 Zr, (b) 120 Sn, and (c) 208 Pb. The calculation 
has been performed with the SkM* parameter set, and 7 = 1 MeV [50]. The experimental data (symbols) are 
taken from Refs. [73-75]. 



For heavier nuclei, the calculation better agrees with experiments [50]. Calculated pho- 
toabsorption cross sections in spherical nuclei 90 Zr, 120 Sn, and 208 Pb are compared with 
experimental data in Fig. 6. The calculated GDR peak shows a splitting, however, this may 
be due to the spurious effect coming from the box discretization. Except for this splitting, 
the results agree well with the experimental data. A single Lorentzian fit for the photoab- 
sorption cross section gives the GDR peak energies of 16.4, 15.2, and 13.3 MeV for 90 Zr, 
120 Sn, and 208 Pb, respectively. The corresponding experimental values are 16.7, 15.4, 13.6 
MeV, respectively. The GDR peak positions are well reproduced within an error of 400 keV. 
We may conclude that the SkM* functional reproduces peak energies of the El resonances 
in heavy nuclei. 

For heavy nuclei, the spreading width was supposed to be a major part of the total damping 
width [69, 70]. However, the artificial width of 7 = 1 MeV, which is supposed to take account 
of missing spreading effects, reproduces the observed GDR width. Although the total damp- 
ing width is about 4 MeV for these nuclei, the spreading width is less than half of the total 
width. In fact, the fragmentation of the strength into non-collective lp-lh states (Landau 
damping) is significant in the present calculation. Thus, the small spreading width (about 
1 MeV) is able to reproduce a broadening of the experimental strength distribution. This is 
also consistent with other recent calculations [71, 72]. 

5.2. Giant resonances in the superfluid phase 

Inclusion of the pair density for systems with superfluidity is, theoretically, a straightforward 
extension of the TDKS to TDKSB equation (See Sec. 2.7). However, in practice, it costs a 
significant increase in numerical task. For instance, the number of matrix elements in Eq. 
(163) is roughly proportional to M 2 for the normal system and M 4 for superfluid systems, 
where M is a dimension of the single-particle model space. At present, it is difficult to adopt 
the 3D coordinate-space representation in Sec. 5.1 for superfluid nuclei [79]. In this section, 
we adopt the symmetry restrictions on the shape of the potentials to reduce the numerical 
costs. 
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Fig. 7 Photoabsorption cross sections in (a) Nd and (b) Sm isotopes as functions of 
photon energy [77]. The Skyrme energy functional with the SkM* parameter set was used. 
The experimental data [80, 81] are denoted by filled squares. 




In the linear response equation (77), the system can oscillate without the external pertur- 
bation v\(oS) = (Vi(u;) = 0), at the eigenfrequencies uj = uj n , Thus, the normal modes of 
excitation are obtained by solving the eigenvalue equation: 

/ vn \ 

(163) 

Solving this eigenvalue equation in the qp basis is the most common method for TDDFT to 
study elementary modes of excitation in superfluid nuclei. Because of numerical difficulties, 
most of current studies with Skyrme EDFs for deformed nuclei are restricted to axially 
deformed nuclei [41, 76-78]. 

We first calculate the qp states in the ground state in the 2D coordinate space assuming the 
axial symmetry. Then, all the necessary quantities are expressed in the qp representation. 
For deformed systems, the number of 2qp states becomes huge and we often need a further 
truncation of the 2qp space [77]. In addition, the residual kernel associated with the long- 
range Coulomb part is neglected in the present calculation. The solutions of the eigenvalue 
equation (163) are obtained using the symmetrization procedure [19]. 

Heavy nuclei with open-shell configurations are supposed to have a superfluid character 
caused by the neutron-neutron and proton-proton pairing correlations. Here, we show results 
for Nd and Sm isotopes. The protons have finite pair densities for all these isotopes, while the 
neutron pair density vanishes for 142 Nd and 144 Sm which correspond to the neutron magic 
number N = 82. These isotopes exhibit typical examples of the quantum phase transitions in 
their ground states, from spherical to prolate shapes, and simultaneously, from the normal to 
superfluid phases, as increasing the neutron number. Actually, the calculated ground states 
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show a spherical shape for 142 ' 144 Nd and ' Sm, and a prolate shape for the others. The 
magnitude of deformation increases as the neutron number changes from 86 to 92. 
In Fig. 7, the calculated photoabsorption cross sections for Nd and Sm isotopes are shown 
together with the available experimental data [80, 81]. The dipole strength at discrete 
eigenenergies are smeared by the Lorentzian with a width 7 = 2 MeV. The GDR peak ener- 
gies well agree with experimental values, and produces the deformation splitting in 150 > 152 Nd 
and 152 ' 154 Sm. The GDR width calculated with 7 = 2 MeV is also in good accordance with 
the experimental values. The nice agreement on the broadening indicates that the smearing 
width 7 = 2 MeV has a good correspondence with the spreading width in these nuclei. 
The isotopic dependence of the peak broadening is well reproduced, even for the transitional 
nuclei. The width for N = 82 and 84 is calculated as T ~ 4.5 MeV, and it gradually increases 
to about 6 MeV for N = 88 ( 148 Nd and 150 Sm), then the peak splitting becomes visible for 
N > 90 and 92. Here, the width T is evaluated by fitting the calculated cross section with a 
Lorentz line. The broadening of the GDR was found to be well correlated with the nuclear 
quadrupole moments [82, 83]. Thus, it is interpreted as the mode- mode coupling effects to 
the low-lying collective modes [70]. In the present calculation, the mode coupling is not 
explicitly taken into account. However, the linear response based on the deformed state may 
implicitly include a part of the coupling effect. Figure 7 shows that the isotopic dependence 
can be well accounted for by the gradual increase of the ground-state deformation. 

6. Large amplitude collective dynamics in shape coexistence 

Low-energy collective modes of excitation in nuclei present unique features of the finite quan- 
tum systems. In contrast to giant resonances discussed in Sec. 5, the linear approximation is 
often insufficient for low-lying collective states in nuclei. The vibrational excitations should 
contain a strong anharmonicity when the stability matrix S has an eigenvalue close to zero, 
namely when the system is close to the critical point of the stability. This kind of situation 
occurs in many nuclei, especially for the quadrupole modes of excitation in a transitional 
situation such as shape phase transition (Fig. 7) and shape coexistence phenomena [84, 85]. 
Therefore, we need to go beyond the linear regime for describing these nuclear phenomena 
of a large amplitude nature. 

One of the unique features of the low-lying collective motion in nuclei is the fact that its 
character significantly changes from nucleus to nucleus. In addition, its structure is affected 
by an interplay between the pairing and deformation correlations, which often results in a 
spontaneous breaking of symmetry. Therefore, it is difficult to introduce an a priori assump- 
tion on the nature of the low-lying collective motion. Therefore, the theory presented in 
Sec. 4 is suitable to find optimal collective manifold which leads to a collective Hamiltonian. 
In this section, we show applications of the theory in Sec. 4. 

In this section, we apply the method in Sec. 4.5 to description of low-lying spectra in neutron- 
deficient 68 Se. This nucleus shows a feature of the shape coexistence: The experimental 
data indicate that there exist rotational bands with different characters, which has been 
interpreted to be those with prolate and oblate shapes [86, 87]. The mixing property of the 
two bands are of significant interest, which influences the excitation spectra and transition 
probabilities. In this section, we show results only for 68 Se. The same analysis on 70,72 Se can 
be found in Ref. [48]. 
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6.0.1. Pairing-plus- quadrupole model. The pairing-plus-quadrupole (P+Q) model is one 
of the most successful models that allows us to describe nuclear phenomena involving the 
quadrupole and pairing degrees of freedom. Baranger and Kumar studied the quadrupole 
motion in the P+Q model, assuming that the collective coordinates are given by the 
quadrupole deformations 7), and that the collective mass parameters are given by the 
cranking formula with phenomenological corrections [88]. However, a study of the same 
model [44] reveals that, even at the minimum point of the potential, the self-consistent 
mass parameters and the property of normal modes are very different from those utilized by 
Baranger and Kumar. In the followings, we also adopt the Hamiltonian similar to the P+Q 
model and study the large amplitude collective motion in Se nuclei. 
The P+Q model in the present study includes the quadrupole pairing and is given by 



k r=n,p 

-\ E ^E^S + ^SV^E^^, 



(164) 



2 

T=n,p K K 



where e& are the spherical single-particle energies and the index k denotes the set of quantum 
numbers (njlm). The monopole pairing operator A^ T ^ , the quadrupole pairing operator 
and the mass quadrupole operator £)' are defined by 

i«t= E 44> B$ = J2 D 2K (M)clcl 

{k >o)er Kier (i65) 

d 2 k =E- p 2A-(fco4 c '' 

k,l 

where D2x(kl) are the (dimensionless) quadrupole matrix elements, modified according to 
a prescription given by Baranger and Kumar [88] . The model parameters are determined by 
adjusting the potential energy surface calculated with the Skyrme energy functional. Those 
values are found in Ref. [48]. 

6.0.2. Results of MFHE. The Hartree-Bogoliubov (HB) approximation to the present 
P+Q Hamiltonian leads to two minima (HB states) in the potential energy surface, in 68 Se: 
The lowest minimum has an oblate shape (/3 = 0.3) and the other minimum with a prolate 
shape (/? = 0.26) lies at energy about 400 keV higher than the oblate one. The obtained 
collective coordinate q 1 has a good one-to-one correspondence to the triaxial deformation 
parameter 7 in the quadrupole deformations (/3, 7) defined by 

/3cos 7 = (ttfo 1 )!^!*^ 1 )), /3sin 7 = (^(q r )\(D 22 + ^^(q 1 )) / ^2. (166) 

This is seen in the top-right panel in Fig. 8. In Fig. 8, a variety of quantities calculated on 
the collective path are plotted as functions of 7. The mass parameter is defined with respect 
to the geometrical length s in the (/?, 7) plane, ds 2 = \J d/3 2 + f3 2 d^ 2 . 

-1/ i\ ( ds V oil _ J ( d P\ 2 , 02 ( rf 7 X 



In the prolate and oblate minima, the lowest modes of excitation in the MFHE correspond to 
the gamma vibrations, and the second lowest to the beta vibrations. The frequency ui, then, 
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Fig. 8 Calculated quantities on the collective path for b8 Se as functions of the triaxiality 7 
[48]. The axial deformation /3(q), the canonical collective coordinate q, the monopole pairing 
gaps A^\q), the quadrupole pairing gaps A^(q) and A^iq), the collective potential V(q), 
the collective mass M(s(q)) } the lowest two eigenfrequencies of the MFHE u) 2 (q), and the 
rotational moments of inertia Ji{q). 



turns into imaginary in the triaxial region (10° < 7 < 50°). The imaginary frequency causes 
no problem in the solution of the MFHE. We also calculate the moments of inertia by solving 
the Thouless-Valatin equations [89] on the collective path, which gives Jk{q), k = x,y,z. 
These are also shown in the bottom-right panel in Fig. 8. They have a characteristic feature 
similar to the moments of inertia of the irrotational fluid [19]. 

6.0.3. Collective Hamiltonian. In the total kinetic energy, the position-dependent rota- 
tional energy is added to the one of the one-dimensional shape vibration described by the 
coordinate q 1 . 

T = ^B i i(q 1 )p iPj (168) 

where i, j = 1, • • • ,4. Here, (p2,P3,P4) are the total angular momentum (I x , I y , I z ), and B ld = 
S l j(B u , JZjT 1 , Jy X 1 Jz 1 )- The kinetic energy term is requantized by means of the Pauli 
prescription: 



2^' ^ " dq v v 71 v " dqi 

1 oii fn u d 1 dG d v _g_ 

2 8q 1 Kq 'dq l Adq 1 Gfg 1 ) 8q l ^ 2J k {q 1 Y K ' 

k=x,y,z 
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where G(q) = Bn{q 1 )J x {q 1 )J y {q 1 )J z {q l ) is the determinant of the metric Bij(q l ). The three 
components Ik of the angular momentum operator are defined with respect to the principal 
axes (x, y, z) associated with the intrinsic (moving-frame) state ^(q 1 )). 
The collective Schrodinger equation is thus given, with q 1 replaced by q hereafter, as 

(T + V{q))V IM n(q, fi) = E Itn V IMn (q, 0), (170) 

where f/Mn(9, ^) represents the collective wave function in the laboratory frame. It is a 
function of the collective coordinate q and the three Euler angles Q, and specified by the 
total angular momentum /, its projection M on the laboratory z-axis, and the index n 
distinguishing different quantum states having the same / and M. Using the rotational wave 
functions T> I MK (^), the collective wave functions in the laboratory frame is written as 

I 

*lMn(q, ty = J2 3>lKn(q)(n\IMK), (171) 
K=0 



(Q\IMK) =J ^J^+s^ VmkW + (-)X^(^)) (172) 

where the sum in Eq. (171) is restricted to even K. Here, <&iKn{q) represents the shape 
vibrational motion described by the coordinate q. 

Normalization of the vibrational part of the collective wave functions is given by 

I 

dr' £ ** IKn (Q)*iKn> (q) = <W (173) 

K=0 

where the volume element is 



dr = dr'dQ = y/\G(q)\dqdn. (174) 

The boundary conditions for the collective Schrodinger equation (170) can be specified by 
projecting the obtained collective path onto the (/?, 7) plane and by using the symmetry 
properties of the Bohr-Mottelson collective Hamiltonian [91]. 

6.0.4- Discussion: Rotational localization. In the left panel of Fig. 9, excitation spectra 
and B{E2) values calculated for 68 Se are displayed together with experimental data. The 
calculation yields two rotational bands, which qualitatively agrees with experiment. The cal- 
culated ground and excited bands have oblate and prolate characters, respectively, however, 
there is strong shape mixing, especially at low spins. Note that the calculated 0^~ state is 
located above the 2^ state. We find that this significant rise of the 0^ excitation energy is 
due to the oblate-prolate shape mixing in the + states, which is much stronger than that 
of the 2 + states. It would be very interesting to identify the 0^" experimentally. 8 
The strong shape mixing in the + states can be confirmed by examining the vibrational 
wave functions displayed in the right panel of Fig. 9. Since the potential barrier between the 
two minima is only 400 keV high along the collective path, this strong mixing is reasonable. 
In fact, the unusual behavior of the excited + state suggests an intermediate situation 



8 In recent calculation with a (2+3)-dimensional collective Hamiltonian, this state is predicted as 
the third + state. 



42/47 



10 



68 



Se 



ASCC 






0.75 




0.5 








0.25 









0.5 


EXP 


0.25 









0.5 




0.25 


6 + 







0.5 


4 + 




2 + 


0.25 








ot stale 


O2 state 




2f stale 


2t state 




4i state 


4i state 


6f state 


K=0 

6J state K-2 

K=A - 

_ K=6 

SUM 



0.75 
0.5 
0.25 


0.5 
0.25 



0.25 




10 20 30 40 50 10 20 30 40 50 60 



Fig. 9 Calculated spectra and B(E2) J, in units of e 2 fm 4 (left) and the collective wave 
functions as functions of 7 (right). The figures are taken from Ref. [48] 



between the oblate-prolate shape coexistence and the 7-unstable model by Wilet and Jean 
[92]. 

It is quite interesting to notice that the shape mixing becomes weak as the angular momen- 
tum increases. The collective wave functions of the the 4 + and 6 + states tend to localize in the 
region near either the oblate or the prolate shape. Namely, it becomes more appropriate to 
characterize the 4 + and 6 + states as oblate-like or prolate-like. We have analyzed dynamical 
origin of this trend and found that the rotational energy plays a crucial role in determining 
localization of the collective wave function. Therefore, this effect may be called "rotational 
localization of collective wave function" or "rotational hindrance of shape mixing" [48]. 

6.1. Further developments 

6.1.1. Multi- dimensional collective submanifold. The ID collective path for 68 Se is 
obtained by following the lowest eigensolution of the MFHE. Figure 10 shows an embedded 
collective path in the (/?, 7) plane. However, the frequency of the second lowest solution is 
only 1 ~ 1.5 MeV higher than the lowest one (Fig. 8). Thus, the extension from ID to 2D 
may be important. We often encounter similar situations which suggest importance of the 
multi-dimensional collective submanifold. Numerically, this is a challenging task, because we 
need to search for self-consistent solutions of the MFHE in the 2D or higher-dimensional 
hypersurface. The technical developments for the multi-dimensional collective submanifold 
is an important future subject. 

In recent papers [93-95], the collective Hamiltonian with the 2D shape degrees has been con- 
structed assuming the one-to-one correspondence between the collective coordinates (q 1 , q 2 ) 
and (f3, 7), combined with following approximations: 

(1) The collective submanifold is determined by the minimization with respect to the mass 
quadrupole operators, neglecting the self-consistency between Eqs. (155) and (156). 

(2) In the MFHE, the curvature terms are neglected. 
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The eigensolutions of the MFHE are used to calculate the collective mass parameters, B 13 = 
f* a f J pB a P. The inclusion of the 2D shape degrees of freedom turns out to be qualitatively 
consistent with the ID calculation in Fig. 9, but further improve the results in comparison 
with experiments [93]. 

6.1.2. Applications to modern energy Junctionals. 
It is also highly desired to apply the method to 
realistic modern energy functionals, such as Skyrme, 
Gogny, and covariant energy functionals. This requires 
a developments of the code for MFHE without the 
axial symmetry restriction. 

Adopting the similar approximations mentioned 
above, there are numerous recent works, in the Skyrme 
energy functionals [96, 97] the Gogny energy func- 
tionals [98, 99], and the covariant EDFs [100, 101], 
to construct the collective Bohr Hamiltonian with 
the microscopic inputs. However, the collective mass 
parameters are further approximated by the Inglis- 
cranking formula. Since this cranking approximation 
has a well-known defect leading to a disagreement 
with experiments, the phenomenological scaling to 
increase the collective mass by roughly 30 — 40 %, 
is often adopted. This comes from the fact that the 
cranking formula neglects the time-odd mean fields 
[19] . Therefore, the replacement of the cranking mass 
by the MFHE mass should solve this problem. 
The MFHE with the Skyrme energy functional has been applied to studies of Cr isotopes 
with the above approximations and a restriction to the ID shape degree [102]. In a very 
recent work [103], a hybrid model, combining the covariant EDF method with the MFHE in 
the P+Q model, has been investigated. Namely, the potential energy surface is calculated by 
the covariant EDF and the collective mass parameters are estimated by the help of MFHE 
in the P+Q model. The model shows a significant improvement in the low-lying spectra in 
7-soft nuclei of Xe and Ba isotopes. These results suggest a promising future of the MFHE 
in the energy functional approaches. 

7. Summary 

Intensive studies in density functional theories (DFT) in recent years have produced numer- 
ous new results and new insights into nuclear structure. It is beyond the scope of the present 
paper to review all of these developments. In this paper, we mainly focus our discussion 
on the basic concepts of the nuclear DFT and applications of the time-dependent DFT 
(TDDFT). 

The nucleus is a self-bound isolated system without an external potential. This produces a 
situation different from many-electron systems with the external Coulomb potentials. There- 
fore, we need modify the original arguments of the DFT, such as the Hohenberg-Kohn 
theorem. In Sec. 2, we presented a justification based on the Hohenberg-Kohn theorem 




Fig. 10 The collective path 
for 68 Se in the (/3, 7) plane. The 
contour plot shows the result 
of the energy minimization with 
the constraints on the mass 
quadrupole operators, D20 and 

D22 + D2-2- 
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modified for the wave packet. It can be formulated with the Kohn-Sham scheme as well. 
To incorporate effects of the nucleonic pair condensation, the density functional should be 
generalized to include the pair density (abnormal density) in addition to the normal den- 
sity. This can be done with the generalized scheme by Bogoliubov. The TDDFT was also 
established, in the same manner as the Hohenberg-Kohn, by the one-to-one correspondence 
between the potential and density depending on time. 

Most applications of the TDDFT are performed in the perturbative linear response regime. 
Even in the linear regime, because of the complexity of the EDFs of nuclei, it is computation- 
ally very demanding to perform the numerical calculations and requires significant works for 
its program coding. There are several approaches to the linear-response calculations, each of 
which has an advantage and disadvantage. This has been discussed in details in Sec. 3. As 
a feasible approach, we have proposed the finite amplitude method (FAM). The FAM only 
requires a minor modification of the existing program of the stationary calculation for the 
nuclear ground-state properties. 

The nucleus is also known to show many collective phenomena at low energies which are not 
able to be described by the linear-response theory. A typical example is given by the fission of 
heavy nuclei, and by the shape phase transition and the shape coexistence. For this purpose, 
we have presented a theory of a decoupled collective space inside the large TDDFT phase 
space in Sec. 4. In case that the collective motion of interest is slow relative to the other 
intrinsic degrees of freedom, self-consistent solutions of the moving-frame harmonic equation 
(MFHE) with the constraint minimization of the EDF provide microscopic determination of 
the collective variables and the collective Hamiltonian. 

In Sec. 5 and Sec. 6, we presented selected results of recent studies with the TDDFT on 
nuclear dynamics in the linear-response regime and beyond the linear regime. Properties 
of giant resonances in nuclei, especially those of the giant dipole resonance (GDR), have 
been studied with the Skyrme energy functionals. Currently available EDFs are able to 
reproduce the GDR in heavy nuclei, however, have a problem for that in light nuclei. We 
also presented the applications of the MFHE to shape coexistence phenomena in 68 Se. The 
calculation suggests a rotation-induced localization of the collective wave functions, which 
produces the oblate-prolate shape coexistence in this nucleus. 

An interplay between theory and experiment has been and will be providing deeper under- 
standing of the nuclear quantum many-body problem. To develop a comprehensive predictive 
theory of the nucleus to answer fundamental scientific questions, still significant components 
are missing from our current understanding. Some of these missing parts can be addressed 
only by the large-scale numerical simulations of nuclear systems. The petascale- and of future 
exascle-computing platforms is expected to pave the way to this goal. 
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